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ABSTRACT 



Aims. We present a new numerical code, ECHO, based on an Eulerian Conservative High Order scheme for time dependent three- 
dimensional general relativistic magnetohydrodynamics (GRMHD) and magnetodynamics (GRMD). ECHO is aimed at providing a 
shock-capturing conservative method able to work at an arbitrary level of formal accuracy (for smooth flows), where the other existing 
GRMHD and GRMD schemes yield an overall second order at most. Moreover, our goal is to present a general framework, based 
on the 3+1 Eulerian formalism, allowing for different sets of equations, different algorithms, and working in a generic space-time 
metric, so that ECHO may be easily coupled to any solver for Einstein's equations. 

Methods. Our finite difference conservative scheme previously developed for special relativistic hydrodynamics and MHD is here 
extended to the general relativistic case. Various high order reconstruction methods are implemented and a two-wave approximate 
Riemann solver is used. The induction equation is treated by adopting the Upwind Constrained Transport (UCT) procedures, appropri- 
ate to preserve the divergence-free condition of the magnetic field in shock-capturing methods. The limiting case of magnetodynamics 
(also known as force-free degenerate electrodynamics) is implemented by simply replacing the fluid velocity with the electromagnetic 
drift velocity and by neglecting the matter contribution to the stress tensor. 

Results. ECHO is particularly accurate, efficient, versatile, and robust. It has been tested against several astrophysical applications, 
like magnetized accretion onto black holes and constant angular momentum thick disks threaded by toroidal fields. A novel test on 
the propagation of large amplitude circularly polarized Alfven waves is proposed and this allows us to prove the spatial and temporal 
high order properties of ECHO very accurately. In particular, we show that reconstruction based on a Monotonicity Preserving filter 
applied to a fixed 5-point stencil gives highly accurate results for smooth solutions, both in flat and curved metric (up to the nominal 
fifth order), while at the same time providing sharp profiles in tests involving discontinuities. 

Key words. Plasmas - Magnetohydrodynamics (MHD) - Gravitation - Relativity - Shock waves - Methods: numerical 



1. Introduction 

Compact objects like black holes and neutron stars interacting 
with the relativistic plasma in the surrounding regions are be- 
lieved to be responsible for many of high energy phenomena in 
astrophysics. The most luminous sources, namely active galac- 
tic nuclei or gamma-ray bursts, are likely to be powered by 
the conversion of gravitational energy of rotating black holes 
into electromagnetic fields and a plasma of relativistic particles 
(Bla ndford & Znaiekll977l) . A similar mechanism had been pre- 
viously proposed to generate the magnetospheric plasma and 
ultimately a Poynting flux domin ated wind from rotating neu- 
tron stars (iGoldreich & Julianlfl969h . The presence of the mag- 
netic field is crucial in all the situations outlined above. The 
magnetic field could also be important in the phases of gravi- 
tational collapse that then give rise to the compact objects them- 
selves, because the freeze-in condition valid for highly conduct- 
ing plasmas would allow an initially negligible field to be en- 
hanced by the collapse to such high intensities to be ultimately 
dominant. The physical frameworks in which these mechanisms 
are treated are usually that of general relativistic magnetohydro- 
dynamics (GRMHD) or, when the electromagnetic field contri- 
bution is dominant over the matt er contribution, that of force- 
free degenerate electrodynamics (Komissarov 2002, 2004), also 



known as magnetodynamics (GRMD, Komissar ov et al.l l2006). 
In both cases the electromagnetic fields interact strongly with 
the plasma, in such a way that freely moving charges are sup- 
posed to screen efficiently any local electric field and to maintain 
quasi-neutrality. 

A great impulse to the study of these complex phenom- 
ena has come from numerical simulations, especially in the 
last decade. Since relativistic magnetized flows are often asso- 
ciated with the formation of strong shocks and different kinds 
of discontinuities, it is thanks to the development of conser- 
vative shock-capturing, or Godunov-type, methods that this 
progress has been possible. After the first applic ations to special 
and general relativistic hydrodynamics (e.g . iFont et afl 1 1994 



1997; 



Eulderink& Mellem al Il994t iBanvuls etal 
19991). see alsolMartf & Mulled d2003l) : lFonll ( l2003l) for reviews, 
Komissarov (1999) first proposed a multi-dimensional shock- 
capturing code for special relativistic MHD (RMHD). These 
schemes are all based on the so-called Roe-type methods, widely 
used in computational gas dynamics, in which the solution of 
the local Riemann problem at any cell interface is constructed 
by means of a full decomposition into characteristic waves. 
However, while this approach is perfectly feasible for purely hy- 
drodynamic flows, in RMHD the spectral structure of the system 
is much harder to resolve, due to the increase in number (from 
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five to seven) and complexity (eigenvalues are to be found nu- 
merically) of the characteristic waves, and to the presence of a 
preferential direction that may lead to non-strict hyperbolicity 
of the local system. Furthermore, the solenoidal constraint for 
the magnetic field in multi-dimensions requires a special numer- 
ical treatment, which must be compatible with the conservative 
approach. 

Within the family of shock-capturing conservative 
schemes, a different strategy was followed in our previ- 
ous investigations on numerica l relativistic hydrodynamics 
dDel Zanna & Bucciantinil |20 02), hereafter Paper I, and 
MHD (iDel Zanna et all 120031) . hereafter Paper II, rely- 
ing on the promising resul ts obtained for classical MHD 
dLondrillo & Del Zannall2000l) . As shown in these works, accu- 
rate and robust schemes can be devised even if the characteristic 
spectral decomposition of the equations is not fully known, or 
exploited, because this lack of knowledge is somehow com- 
pensated by resorting to higher (third) order reconstruction of 
intercell variables, leading to a more accurate setup of the local 
Riemann problem. By doing so, even simple one or two-wave 
approximate Riemann solvers (also known as central-type 
schemes) are capable of resolving all kinds of structures, thus 
avoiding the problems related to the complexity in spectral de- 
composition at the price of a slightly higher numerical diffusion 
of contact and Alfvenic discontinuities. Many other shock- 
capturing numerical codes for RHMD and GRMHD (some of 
them even with evolving space-time met ric) share the same 
philosophy of a simplified Riemann solv er (jGammie et al. 2003; 
Duez et al.l 120051; IShibata & Sekiguc hU 120051; iLeismann et alj 
20051; iMignone & Bodol 120061; lAnton et alj 120061) . though all 



of them are based on finite di f ferenc e or finite volume second 
order schemes. In lAnton et al.l ((2006) an RMHD Roe solver is 
also used in s ome tests, via a l ocal coordinate transformation 
to flat metric dPons et al.lll998l) . Moreover, different methods 
other than Goduno v -type have also been proposed for GRMHD 
dKoide et all I l999t iKoidd 120031; iDe Villiers & Hawlevl 120031: 
lAnninos et a l.l l2005l) and (GR)MDlSpitkovskv[l2006l) See also 
the reviews bv lWilson & Mathews! d2003l) : lFontl(l2003l) . 

These codes have been extensively applied to many as- 
trophysical situations involving relativistic plasmas and com- 
pact objects. Relevant examples of these applications in- 
clude the validation of the Blandford-Znajek mechanism 
for th e extraction of ro ta tional e nergy from a Kerr black 
hole dKomissarovl I200l1 iKoidd l2003t iKomissarovl 120041: 



McKinnev & Gammie 2004; Komissarov | 12005b iMcKirmevI 
2005); the spin evolution of a black hole under the ef- 



fect of different physical processes (Gam mie et al.l [2004); the 
probl em of jet formation in a black hole-accre t ion disk sys- 
tem dKoide_eTaII 120001: foe Villiers et ail l2003j iMizuno et al.l 



2004 INishikawa et alj|2005t IDe Villiers et al.ll2005t IMcKinnevI 



2006bj lHawlev & KroliM l2006t iKoide et all 120061) ; the time 
evolution of a neutron star m agnetosphere, both in the MHD 
regime (Komissarov 2006b) and in the forc e-free approxi- 
mation (McKinnev! l2006ct ISpitkpvskvl 20061); the acce lera- 
tion of magnetized pulsar winds (Bucciantini et al. 2006) and 
the d ynamics and emission properties of their related neb- 
ulae (Komissarov & Lyubarsky 2004; |Del Zanna et al.l [20041; 
iBucciantini et al.ll2005t IDel Zann a et al.ll2006l) : the morphology 
and the dynamics of axisymm etric relativistic jets w ith differ- 
ent magnetic field topologies (Leismann et al. 1 (20051) : the col- 
lapse, in full genera l relativity of a hyp er-massive neutron star 
( Shib ata et al.l 120061: iDuez et aLll2006al). also including the ef- 
fects of differential rotation (IDuez et al.ll2006bl) . All of these ap- 
plications, that do not pretend to provide a complete list, surely 



give a sample of the fundamental contributions that numeri- 
cal simulations have been offering to our understanding of the 
highly complex physical processes induced by the relativistic 
plasma around compact objects. 

In this paper we present the main features of our new 
GRMHD code ECHO, based on an Eulerian Conservative High 
Order scheme, that completes and refines our previous works for 
special relativity (Paper I and II). The issue of high numerical 
accuracy in conservative schemes becomes of great importance 
when not only shocks and discontinuities, but also fine smooth 
structures like turbulent fields and waves, are of primary interest. 
These small scale structures can be smeared out by the excessive 
numerical diffusion typical of low order schemes. Furthermore, 
higher than second order accuracy is desirable when moving 
to 3-D, where numerical grids are necessarily limited in size. 
This specially applies to GR, due to the gradients of the metric 
terms that must be treated with appropriate resolution. High or- 
der schemes are commonly used in classical gas dynamics (e.g. 
IShulll997l) . and the general recipes to app l y thes e methods to 
MHD were given in lLondrillo & Del Zannal ( I2000ll2004l) . where 
the solenoidal constraint for the magnetic field was enforced 
as a built-in condition (Upwind Constrained Transport method, 
UCT). Here we extend this framework to GRMHD by taking ad- 
vantage of the formalism for the 3 + 1 splitting of space-time (e.g. 
iThorne & MacDonaldll9 82). Specifically, we write all terms en- 
tering the conservative form of the GRMHD equations as quanti- 
ties measured by the so-called Eulerian observer associated with 
the three-dimensional metric (not necessarily diagonal), high- 
lighting the closest possible comparison with the equations of 
MHD and RMHD by using three-dimensional vectors and ten- 
sors alone. As a consequence, we are able to write the source 
terms in such a way that they do not contain four-dimensional 
Christoffel symbols explicitly, and are therefore very easy to im- 
plement numerically. We then inc orporate in the 3 + 1 formal- 
ism the modifications proposed by McKinnev (2006a) to allow 
a GRMHD code to solve the equations in the force-free limit of 
magnetodynamics (GRMD). 

The plan of the paper is as follows. In Sect. [2] we present 
the 3 + 1 form of the GRMHD equations. Sect. [3] contains a 
description of the essential features of our numerical scheme. 
Sects. [4]and[5]are devoted to a presentation of the most important 
numerical tests performed in GRMHD and GRMD, respectively. 
Finally, the conclusions are reported in Sect. [6] In the following 
we will assume a signature {-, +, +, +} for the space-time metric 
and we will use Greek letters p, v, A, . . . (running from to 3) 
for four-dimensional space-time tensor components, while Latin 
letters i, j, k, . . . (running from 1 to 3) will be employed for three- 
dimensional spatial tensor components. Moreover, we set c = 
G = Mq = 1 and make use of the Lorentz-Heaviside notation for 
the electromagnetic quantities, thus all V47r factors disappear. 

2. GRMHD equations in 3 + 1 conservative form 

2.1. Covariant approach 

We start with a brief presentation of the GRMHD equations in 
covariant form. Standard derivations of the laws of fluid dy- 
namics and elect rodynamics in covariant f o rm may b e foun d 
in books such as | Landau & L ifshitz (1962); Weinberg! d 19721) : 
Misn er etail d 1973b. while for the MHD equat i ons an d their ba- 
sic properties see iLichnerowiczl dl967l) : lAnild (Il989l) . Consider 
an ideal fluid interacting with an electromagnetic field. The cor- 
responding Euler equations are 

V>m") = 0, (1) 
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V^ v = 0, 



(2) 2.2. The 3 + 1 splitting of space-time 



where V^, is the space-time covariant derivative. Eq. ([TJ is the 
usual mass conservation law, in which p is the mass density as 
measured in the (Lagrangian) frame comoving with the fluid 
four-velocity u 11 . Eq. is the law of momentum-energy con- 
servation, where the total momentum-energy tensor is made up 
by two contributions, T 1 ™ = + T^ v , one due to matter 



(UV 



T% = ph u M u + pg 
and the other due to the electromagnetic field 



= F» A F vA ■ 



(F A "F AK )g> 



(3) 



(4) 



In the above expressions g^ v is the space-time metric tensor, h = 
1 + e + pip is the specific enthalpy (including rest mass energy 
contribution), e is the specific internal energy, p is the thermal 
pressure, F^ v is the electromagnetic field (antisymmetric) tensor. 
When considered separately, the two components of the stress 
tensor are not conserved 



V T /JV - -V - - 1 F^ v 

v n 1 m — V f ~ M ' 



(5) 



where J M is the four-vector of current density and the last term 
is the electromagnetic force acting on the conducting fluid. The 
fields obey the two Maxwell equations 



v^" = -r 



(6) 
(7) 



where F*^ = he^^F^ is the dual of the electromagnetic ten- 
sor, and €^ vAk is the space-time Levi-Civita tensor density, that 
is ei"** = (-g)- l/2 [pvAK] (and e, vAK = -(-g) l ' 2 [pvA«]), with 
g = det{gjj V } and [pvAi<] is the alternating Levi-Civita symbol. 

Since we are dealing with a (perfectly) conducting fluid, a 
general relativistic extension of (ideal) Ohm's law is needed. 
This translates in a condition of vanishing electric field in the 
comoving frame 



F" v u v = 0. 



(8) 



From a physical point of view it means that the freely moving 
charges in a plasma are supposed to be always able to screen any 
electric field that may arise locally. The extra condition imposed 
on F 1 " in Eq. dS) makes the first Maxwell equation redundant, 
and Eq. © is only needed to calculate the four-current J 1 *, which 
is now a derived quantity like in non-relativistic MHD. The sys- 
tem of GRMHD equations is then closed by choosing an equa- 
tion of state (EoS) p = p(p, e). Different relativistic EoS may be 
employed, and thus we will leave it unspecified in our formula- 
tion. However, all numerical tests presented here will make use 
of the standard y-law for a perfect gas 



p(p,e) = (y- l)pe=> h = 1 + 



r p 



y- 1 p 



(9) 



with y — 5/3 for a non-relativistic fluid and y — 4/3 when p » p 
(ph — > Ap). Finally, note that for an ideal fluid (thus in the ab- 
sence of shocks or other sources of dissipation) the total energy 
conservation law is equivalent to the adiabatic equation 



M *%S = : 



Vfiipsu' 1 ) = 0, 



(10) 



even in the GRMHD case (e.g. lAniielH989h . Here s is any func- 
tion of the specific entropy (in the comoving frame), and in the 
case of a fluid with a y-law EoS we can take s = pip 7 . 



In spite of their elegant and compact form, the GRMHD co- 
variant equations described above are not suitable for numeri- 
cal integration, where the temporal coordinate must be clearly 
singled out. The most widely used formalism is that based on 
the so-called 3 + 1 decomposition of the equations. For a com- 
prehensive treatment a nd ref erences the reader is referred to 
Thorne & MacDonald (1982), or, for a more recent work, see 
iBaumgarte & Shapirol y003). 

In the 3 + 1 formalism, the four-dimensional space-time is fo- 
liated into non-intersecting space-like hyper-surfaces "£,, defined 
as iso-surfaces of a scalar time function t. Let then 



(n M n" = -l) 



(ID 



be the future-pointing time-like unit vector normal to the slices 
"L,, where a is called the lapse function. The observer movin g 
with four- velocity n A< is called Eulerian (ISmarr& Yorkl [l978>. 
and all quantities may be decomposed in the corresponding 
frame. Thus, any vector V (or similarly a tensor) may be pro- 
jected in its temporal component V" = —n fl V fl and spatial com- 
ponent ± V — (gy +n fI n v )V v . In particular, a three-dimensional 
spatial metric y^ v can be induced on S, by the four-dimensional 
metric. Application of the projection operator gives 



TV =J - 8ltV = 8fIV + Ifttly, 



(12) 



so that we can also identify ±=±y= y^. At this point, it is con- 
venient to introduce a coordinate system x M = (t, x') adapted to 
the fol iation The line elem ent is usually given in the so-called 
ADM dArnowittet al.ll 1 9621) form: 



ds z =-a z dt z +y ij (dx' + /3'dt)(dx J + /3 J dt), 



(13) 



where fi^ is called shift vector, an arbitrary spatial vector 
(/S^Hju = 0). Notice that the spatial metric y (/ can now be used for 
the raising and lowering of indices for purely spatial vectors and 
tensors. In this coordinate system the unit vector components are 



(-«,0,), 



{\ I a,- 0/a), 



(14) 



and any spatial vector V (or tensor) must necessarily have a 
vanishing contravariant temporal component V' = 0, whereas its 
covariant temporal component is V, = g^tV = y6,V", in general 
different from zero. The gradient of the unit vector can also 
be split into spatial and temporal components as follows 



-K, 



(15) 



where is the extrinsic curvature of the metric (a spatial sym- 
metric tensor) and a v is the acceleration of the Eulerian observer 
(a sp atial vector too). Finally, it is possible to demonstrate that 
re.g. lYorklll979h 



V v In a, 



(16) 



another property that will be used later on. 

The next step is then to decompose all quantities appearing 
in the GRMHD equations of Sect. 12. ll into their spatial and tem- 
poral components. Hence, we define 



u" = r«" + rv", 

= W y + S"n v + n> 1 S v + Un"n v , 
pnv _ n v E v _ &n v + e ^« Bx n K , 

pit* _ n M B v _ B n n v _ e ^ Exn ^ 



(17) 
(18) 
(19) 
(20) 



L. Del Zanna et al.: ECHO: an Eulerian Conservative High Order scheme for GRMHD and GRMD 



4 

where all the new vectors and tensors are now spatial and corre- 
spond to the familiar three-dimensional quantities as measured 
by the Eulerian observer. In particular is the usual fluid ve- 
locity vector of Lorentz factor F = u , for which 

V = ii'/T + 0/a, r = a M ' = (l-v 2 r 1/2 , (21) 

where v 2 = V;V' and we have used the property u^ = — 1. An 
alternative quantity, u'/u' = av' -B\ usually referred to as trans- 
port ve locity, is sometimes used in stead of the Eulerian velocity 
v' (see Baum garte & Shapiroll2003l). The definition in E g. d2"TT l 
agrees with the treatments by | Thprne & M acDonald (1982); 
ISloan & Smarjdl985h :IZhang ( 1989) a nd it is the m ost appropri- 
ate for numerical integration dBanvuls et al. 1997), since in the 
3 + 1 formalism v' is a real three-dimensional vector while u'/u' 
is not. The decomposition of the momentum-energy stress ten- 
sor gives the quantities U = T m , S" =X T^, and F =J- T» y , 
which are respectively the energy density, the momentum den- 
sity and the spatial stress tensor of the plasma. Finally, the spatial 
electromagnetic vectors in Eqs. ( fT9ll20l i are defined as E 1 * = F w 
and W = F*"^, that is, in components 

E' = aF'\ B' = aF* H . (22) 

2.3. Derivation of the 3 + 1 GRMHD equations 

The set of GRMHD equations in 3 + 1 form is derived from 
that in Sect. 12. 1 1 by applying the space-time decompositions of 
Eqs. ( 11711201 ). Here we are interested in retaining the conserva- 
tive form, as needed by any shock-capturing scheme {Fontl2 003 ; 
IShibata & Sekiguchil2005l;lDuez et aLll2005UAnt6n et al.ll2006i 
In this respect, we improve on these works by making use of 
purely three-dimensional quantities alone, in a way to maintain a 
close relation to classical MHD as much as possible and to sim- 
plify the expression of the source terms. By applying standard 
covariant differentiation relations, the set of GRMHD equations 



becomes 








(-gT m d»[{- 


-g) 1/2 pu»] = 


0, 


(23) 


(-gr u %[(- 


-g) U2 T"j] = 


jT^ djgjxv, 


(24) 


(-gr u %[(- 


-g) 1,2 T» v n v ] 


= r"%« v , 


(25) 


i-gr u %[(- 


-g) 1/2 F w ] = 


= o, 


(26) 


i-gr u %[(- 




:0, 


(27) 



where Eqs. ([T), ©, and (0 have been split into their spatial and 
temporal components and the symmetry properties of T /JV and 
F*i™ have been exploited. Eqs. d21H22b must now be plugged into 
the above equations to yield equations for the three-dimensional 
quantities alone. Moreover, it is easy to verify that the source 
terms on the right hand side are split as 

\T^d jgllv = \W ik d jyik + a^SidjP - Udj In a, (28) 

7^%n v = -KijW ij - S j dj In a, (29) 

where the properties of the extrinsic curvature have been used. 
Notice that only spatial derivatives along j appear in Eq. ( f2~8b . 
so that the corresponding flux is a conserved quantity in the sta- 
tionary case. Finally, it is convenient to introduce the standard 
boldface notation for (spatial) vectors and to define V =± V as 



the three-dimensional covariant derivative operator for the met- 
ric jij (providing the familiar divergence and curl operators), so 



that the final form of the GRMHD equations is then 

y- l/2 d, (y l/2 D) + V • (avD - BD) = 0, (30) 

y- 1/2 d, (y 1/2 S) + V ■ (aW-BS) = (V/J) ■ S - UVa, (31) 

y- l/2 d, (y 1/2 U) + V ■ (aS -fiU) = aK :W - S ■ Vor, (32) 

y~ lj2 d, (y 1/2 B) + V x (aE + J3 x B) = 0, (33) 

V • B = 0, (34) 



where y — det{y, j} is the determinant of the spatial metric (not to 
be confused with the adiabatic index), for which (-g) 1 ^ 2 = ay 1 / 2 . 

Let us analyze the above system in detail. Eq. ( 1301 ) is the con- 
tinuity equation for D = pF, that is the mass density measured by 
the Eulerian observer. The momentum equation, Eq. OTb . con- 
tains the divergence of the tensor W, leading to source terms 
present also in MHD and RMHD when curvilinear coordi- 
nates are used, whereas the last term with the gradient of the 
lapse function becomes the standard gravitational force in the 
Newtonian limit. Eq. d32b is the energy equation, in which the 
extrinsic curvature must be evolved through Einstein's equa- 
tions or, for a stationary space-time, it is provided in terms of 
the covariant derivatives of t he shift vector components (e.g. 
Misn er et al ] |1973UYorklll979h . Here we write 

aK : W = \ W ik p j djy ik + W<9 ; p\ (35) 

where again the symmetry properties of W 1 * have been used. 
Eq. (l33l is the GRMHD extension of the induction equation, 
written in curl form by exploiting usual vector calculus relations. 
Note that the (spatial) three-dimensional Levi-Civita tensor den- 
sity e^ vA = e"^- 1 , for which e ijk = y^ l2 [ijk] and e ijk = y l/2 [ijk], 
is implicitly defined in Eq. d33"T l. Finally, Eq. d34b is the usual 
divergence-free condition. Notice that the above treatment is 
valid in a generic system of curvilinear coordinates, not neces- 
sarily under the assumptions of diagonal spatial metric tensor 
or vanishing expansion factor V ■ /? (e.g. Kerr metric in Boyer- 
Lindquist coordinates). In the absence of gravity, that is when 
a - 1, /? = 0, K - 0, and d,y = 0, the above equations reduce to 
the familiar set of RMHD in curvilinear coordinates. 

The expression for the stress tensor, momentum density, and 
energy density in terms of the fluid and electromagnetic quanti- 



ties are, from Eqs. ( 1 17H20b : 

W = phY 2 vv - EE - BB + [p + \(E 2 + B 2 )]y, (36) 

S = p/iF 2 v + ExB, (37) 

U = phY 2 -p+ \(E 2 + B 2 ), (38) 



where we have indicated with the symbol y the spatial metric 
tensor of components yy. The matter and electromagnetic field 
contributions have been expanded by using Eqs. ([3]|4]i written in 
terms of scalars and the spatial vectors v, E, B alone. In the 3 + 1 
split the Ohm relation for MHD in Eq. (|8]l becomes the usual 
freeze-in condition 

E = -v x B, (39) 

that allows us to close the set of GRMHD equations. Note that all 
the above relations, from Eq. (l36l l to (|39l ), are exactly the same 
as in the special relativistic case (though in Paper II a different 
formalism was employed). Moreover, the non relativistic limit is 
found by letting v 2 « 1, p « p, and £ 2 « £ 2 « p. Thus, by 
simply changing the definition of D, W, S, U and by neglecting 
gravity terms (or reducing them to the Newtonian limit), one has 
the formal setup of a conservative scheme for classical MHD in 
generic curvilinear coordinates. 
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3. The ECHO scheme 



3.1. Discretization and numerical procedures 



The set of conservative GRMHD equations described in 
Sect. I2.3I may be rewritten in a compact way as follows. The 
five scalar fluid equations are 



3/U + diT' = S, 



(40) 



where the conservative variables and the correspondent fluxes in 
the i direction are respectively given by 



(41) 



and the factors y 1 ^ 2 have been included in the definition of these 
new quantities. In the case of a stationary metric, used in the re- 
mainder of this paper for code testing, the source terms become 





D 






= r 1/2 


Sj 


, r = r 1/2 


aW)-f?Sj 




u _ 





M 2 





\aW k djj ik +SidjF - Udja 
\W k ^djj ik + Wiidjft - S J t),a 



(42) 



in which the extrinsic curvature in the energy equation Eq. 
has been replaced by the derivatives of the metric according to 
Eq. d35l l. As far as the induction equation is concerned, it is con- 
venient to introduce the new quantities 



& = y l,2 B<, 

Si = aEi + e ijk B j B k = -[ijkW j S k , 



(43) 
(44) 



where 'V-i - av' - B' is the transport velocity. Eq. (l33l may be 
then rewritten in the form 



80 + [ijk]djS k = 0, 



(45) 



and the related non-evolutionary constraint Eq. ( l34l . expressed 
in terms of the new variables 3', simply becomes 



80 = 0. 



(46) 



Notice that, thanks to our definitions, Eqs. (1401 . d45l l. and (|46*T i 
retain the same form as in Cartesian coordinates (with external 
source terms). Eq. ( |45T > is the conservation law for S', which dif- 
fers from the form of Eq. ( l40l i. basically due to the antisymmetric 
properties of the Faraday and Maxwell tensors. The curl nature 
of the induction equation and the divergence-free constraint must 
be maintained in the numerical scheme by employing consistent 
algorithms. 

In the following we describe the numerical procedures em- 
ployed in our new ECHO code. The scheme is quite general and 
can be applied to any set of physical laws with evolution equa- 
tions in the form of Eqs. ( |40H45l l, with the additional constraint 
of Eq. d46t : physical modules are available for classical MHD, 
special RMHD, GRMHD, and GRMD (see Sect.[3T4l>. The gen- 
eral recipes for the correct treatment of the divergence-free con- 
dition in any shock-capturing MHD-like scheme, regardless of 
the discretization technique (finite volume or finite difference), 
accuracy order, int erpolation methods, and Riema nn solver, have 
been presented in lLondrillo & Del Zannal ([2004). That method 
was named Upwind Constrained Transport (UCT) and here we 
follow its guidelines. In particular we will adopt the same build- 
ing blocks already employed in Paper II, namely finite difference 
discretization, high order component-wise reconstruction meth- 
ods (additional algorithms will be proposed here), a two-wave 
approximate Riemann solver, and multi-stage Runge-Kutta for 
time integration. 



The starting point is the discretization of the GRMHD equations. 
Here we assume a finite difference approach and thus we adopt 
the corresponding version of UCT. This is known to be more 
convenient than finite volume methods for high order treatments 
of multi-dimensional problems , since only 1-D reconst ruction 
algorithms are needed (e.g [Shul[T997tlLiu & Osheril 19981) . Let r 
be the order of spatial accuracy requested for the scheme. Given 
a computational cell of edge sizes hi, the fluid conservative vari- 
ables Kj are defined at cell centers C with a point value repre- 
sentation, that is Uj is the numerical approximation, within an 
accuracy r, of the corresponding analytical function. The other 
conservative variables are the &' components, which are here 
discretized as point values at cell interfaces Sf, normal to di- 
rection i. This discretization technique is k nown as st aggering, 
first introduced for Maxwell's equations by Yee ( 1966 j) and later 
applie d to the GRMHD induction equation by IE vans & Ha wlev 
( 1988). In a conservative approach, the spatial differential opera- 
tors of divergence and curl are translated numerically by making 
use of the Gauss and Stokes theorems, respectively. Fluid fluxes 
T'j are to be calculated at cell faces St, while magnetic fluxes 

& k mu st be calculated at cell edges L t, parallel to the direction 
k (see lLondrillo & Del Z anna 2004). The spatially discretized 
GRMHD equations are then written in the following way 



At 

I' 



m s ; + J][ijk]^([S k ] L; - [8 k ] L -) = 0, 

j,k I 



(47) 



(48) 



known as semi-discrete form, since the time derivatives are left 
analytical. Here the hat indicates high order approximation of the 
numerical flux function, as it will be described at steps 4 and 8 
below, and we have indicated with ± the opposite faces, or edges, 
with respect to the direction of derivation. Time evolution is here 
achieved by means of Runge-Kutta integration schemes. In the 
same framework, the non-evolutionary solenoidal constraint be- 
comes 



(49) 



Given the particular discretization of the conservative quan- 
tities and of their corresponding numerical fluxes, the procedures 
required by the UCT strategy may look rather involved, in par- 
ticular for high order implementations. In the ECHO scheme we 
have made an effort to simplify them as much as possible, espe- 
cially as far as the induction equation and the metric terms are 
concerned. We describe these procedures in the following ten 
steps. 

1. Given the value of the conservative variables at time t, we 
first interpolate the magnetic field components &' from the 
corresponding staggered locations Sf to cell centers C, for 
every direction i. For a second order r = 2 scheme we simply 
use 



[&ic = \(m s - + m st ), 



(50) 



whereas larger stencils are employed for higher order inter- 
polations (see Sect. IA.ll in the appendix). The set of conser- 
vative variables 



W= [t/,S] 3 



(51) 
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is now entirely defined at cell center C. From this we can 
then derive the primitive variables f, that is any set of 
physical quantities such that the functions 1i - t((P) and 
T l = T'CP) are uniquely defined. Here we use 



P=[p,v,p,Bf 



(52) 



for all MHD-like modules in ECHO. In Sect. Owe describe 
the inversion routines implemented for this choice of primi- 
tive variables. 

2. For each direction i, say x, we reconstruct the point value 
approximations of the left (L) and right (R) upwind states of 
primitive variables, from C to S p. 



[P^hi =RfW>j\c\), 



where K 1 ^ is the 



(53) 



1-D reconstruction routine, here named 
REC, applied to a stencil {[!Pj]c) of cell centered values 
along x. The index j runs through all fluid components and 
the transverse magnetic field components. This is because 
the main assumption in UCT is that the longitudinal B x com- 
ponent does not present different upwind states at S At this 
location one can safely assume B xL = B xR = y~ l l 2 !B x . 
In ECHO different reconstruction routines are implemented. 
All of them are treated component-wise, that is avoid- 
ing decomposition into characteristic waves. For schemes 
with overall r — 2 accuracy we may use simple TVD- 
like reconstructions based on limiters (e.g. MM2 for the 
MinMod, MC2 for Monotonized Centered). For r > 2 we 
have a choice of ENO-like routines: ENQ3 f or the third- 
order original ENO method (lHarten et al.ll 19871). C ENQ3 for 
the Convex-ENO scheme by iLiu & Osherl (1998) (see also 
Paper I), WENQ5 f or the Weighted-ENO fifth order scheme 
(lJiang&Shul [l996). Moreover, in the tests of Sect. [4] and 
[5] we will l argely make use of the Monotonicity Preserving 
scheme by Suresh & Huvnhl d 1997b . implemented in ECHO 
as MP5, which is based on interpolation built over a fixed 
5 -point stencil (we recall that adaptive stencils are used in 
ENO schemes), followed by a filter, basically a combina- 
tion of limiters to preserve monotonicity near discontinu- 
ities. Notice that our reconstruction process is based on up- 
wind, non-oscillatory interpolation techniques (thus from 
point values to point values), while in the numerical litera- 
ture reconstruction via the primitive function (or equivalently 
from cell averages to point values) is typically discussed. All 
interpolation coefficients for high order methods are thus dif- 
ferent, and these are calculated in Sect. |A.2l of the appendix. 
The upwind flux for the fluid part is then derived in terms of 
the two-st ate recons tructed primitive variables. In Roe-like 
schemes dRodll98 Tb this task is achieved by a field-by-field 
spectral decomposition of the local Jacobian 7x7 matrix 



J\ x 



dT x 

Ww x 



W x = YU,W,&f 



(54) 



where S x acts like a given parameter in this local 1-D sys- 
tem. The eigenvalues of Sl x , typically calculated at some 
averaged state, provide the speed of each characteristic 
wave. Here we use the HLL approximate Riemann solver 
(lHarten etalJI 19831) which is based on the knowledge of the 
two highest (in absolute value) characteristic waves alone. In 
GRMHD they correspond to the fast magnetosonic waves, 
see Sect. 13.31 If A x ± are the requested speeds, calculated at 
both left and right states, we then define the quantities 



and the HLL upwind fluid flux function is 

a x + T xL + a x _Tf - ala x _(V{* - tlj) 



(56) 



where all quantities are calculated at S\ for each component 

./' and where T xL ' R = T X CP L ' R ), W L ' R = ( U{ r P UR ). At the 
same location we also calculate the upwind transverse trans- 
port velocities and we average them as follows 

— j a x JV jL + a x <Vi R 

V= ~ , j = y,z. 

a, + a 



(57) 



These quantities are saved and will be used at step 6 for the 
calculation of the electric field needed in the induction equa- 
tion. The coefficients a x ± are saved too, since they will be 
needed at step 7 for the magnetic fluxes and at step 10 for 
the timestep definition. Local Lax-Friedrichs is retrieved as 
usual when a x + - a x _. 
4. The numerical fluid flux function is retrieved by means of 
an additional high order procedure, named DER, which al- 
lows one to obtain a high order approximation from the point 
value quantities calculated at the same intercell locations: 



[T*] st = £>A\[Tf ] st }) ■ 



(58) 



This correction step is necessary to preserve the accuracy 
in the calculation of spatial partial derivatives for high or- 
der schemes, while it can be avoided for low order r < 2 
schemes, for which the DER operator is just an identity. In 
the tests with r > 2 presented in Sect. [4] we use fourth or 
sixth order fixed-stencil algorithms (see Sect. |A.3l in the ap- 
pendix). 

5. The fluid flux functions are recovered for all directions i by 
repeating steps 2-4 and the spatial operator in Eq. d47l ) is 
calculated. The source terms [S]c are also worked out so 
that we are ready for the Runge-Kutta time-stepping cycle as 
far as the fluid part is concerned. 

6. The induction equation is treated as follows. Let us concen- 
trate on the magnetic flux [S z ]Lt, the other components are 
found with similar strategies. First we need to reconstruct 
the quantities e V x , < V y , S x , and S y from faces S + and Sy to 
the edge L* t o be combined there in a four- state upwind nu- 
merical flux (lLondrillo & Del Zannal 120041) . Exploiting the 
uniqueness of the numerical representation of as dis- 
cussed at step 2, it is sufficient to reconstruct the following 
quantities 

fV XL '\t = K L x ' R ({VV X h t }), [S vi ' R k+ = K L X - R ({[&'ht }), (59) 



vL,R 



KfiWV ] si }), [S xL ' R ] Lt = ^ S ({[S- C ] S .}), (60) 



where *k (j = x,y) were saved at step 3. 
7. The HLL numerical flux for the magnetic field can be then 
defined as 



S 7 



a x ± = max{0, ±A*(!P L ), ±A X JP R )} 



(55) 



a x + V L S yL + a x JV XR ®> R - alai(&' R - S' L ) 
al + a x _ 

a y JP yL !B xL + a y <V yR !B xR - a y + a y _{S xR - S xL ) 

+ — y-^ > (61) 



which coincides with t he fou r-state formula presented in 
lLondrillo & Del Zannal ([2004). Note that our flux formula 
contains upwinding in the two directions x,y and reduces 
correctly to the expected flux for 1-D cases. 
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8. Following the same strategy as in step 4 the DER operation 
is needed to recover numerical fluxes with appropriate ac- 
curacy. Each magnetic flux component actually requires two 
distinct high order corrections 

[S z ] Lt =£>y({[£ z ]zf}), j = x,y (62) 

as Eq. (08]) contains both x and y differencing of S z , 

9. The spatial derivatives in Eq. d48l are then calculated for 
each direction and also the induction equation is ready for 
time integration. 

10. Runge-Kutta time-stepping can be finally achieved, and the 
whole procedure to update the set of conservative variables 
< W must be repeated for each sub-cycle. Here we use for 
r < 2 the classical Heun (or improved Euler) second or- 
der scheme (RK2), whereas for r > 2 it is convenient to 
use corres pondingly higher o rder methods, like those de- 
scribed in Ishu & Osherf (Il988l) . In ECHO we have imple- 
mented their third order scheme (RK3, see also Paper I). Like 
in all explicit schemes, the timestep Af is limited by the CFL 
(Courant-Friedrichs-Lewy) condition < c < 1 (we will 
always use c = 0.5 in the tests presented) and is defined as 



max,(a^//z,) ' 

where a' M = max({[a' f ]s+}, {[al]s+)) are the maximum 
speeds over the whole domain, for each direction i. Gravity 
contributions to Af are included in the a' M definition via 
the metric terms contained in the GRMHD speeds A' ± (see 
Sect. ED). 

Compared to our previous implementations for classical 
MHD and RMHD, the ECHO scheme presented here is slightly 
simpler. First, the DER operator is now based on fixed, symmet- 
ric stencils, rather than adaptive like in REC (see the appendix). 
As far as the induction equation and the related divergence-free 
constraint are concerned, the use of the magnetic vector potential 
is avoided and the primary magnetic field (staggered) compo- 
nents f or the UCT strategy are now \S']s+, rather than [2?'],?+ 
like in Londr ilio & Del Zannal (|2004), so that magnetic fields 
are also easier to initialize. Moreover, it is easy to verify that 
Eq. d49l is satisfied algebraically at all times regardless of the 
value of r. This is because, when using Eq. d48b in the time 
derivative of the solenoidal condition, the electric field compo- 
nents (now with corrections along the two orthogonal directions) 
cancel each other, due to the commutativity of the DER opera- 
tors applied. Obviously this property holds only for fixed-stencil 
procedures. 

Finally, notice that the metric terms are needed at cell cen- 
ter (where also their derivatives must be given) and at intercells, 
but not at cell edges. This is due to our definitions of the <V 
and £' components, already containing the metric terms needed 
for the calculation of the electric field fij. The components of 
the metric tensor and their derivatives are here provided analyti- 
cally. Another option (e.g. when solving Einstein's equations) is 
to interpolate and derive them, wherever needed, with high order 
procedures as those described in the appendix. 

3.2. Primitive variables 

As we have seen in Sect. [3] in step 1 the primitive variables f 
must be derived from the set of conservative variables 1V at cell 
centers. The problem is exactly the same as in special relativistic 
MHD, that is: 

[D,S,U,B]^[p,v,p,B], (64) 
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with B acting at the same time as a conservative and primi- 
tive variable. Her e we basical l y foll ow the strategy outlined in 
Paper II, see also iNoble et aU d2006l) for further discussion and 
comparison of different techniques. The full system is first re- 
duced to a 2 x 2 set of nonlinear equations in the variables x = v 2 
and y = phT 2 . Let us rewrite Eqs. d37b and d38l using Eq. d39b 
for the electric field, and then calculate S 2 and S ■ B. After some 
simple algebra, the unknown variables may be found by solving 
the system F\ — 0, F2 = 0, where 

Fi(x,y) = (y + B 2 fx - y- 2 (S ■ B) 2 (2y + B 2 )-S 2 , (65) 
F 2 (x,y) = y - p + i(l + x)B 2 - \y~\S ■ B) 2 - U, (66) 

with p = p(x, y) to be specified according to the EoS employed. 
Once x and y are found, the required primitive variables are given 
by the relations 

p = D(l - x) l/2 , (67) 
v = (y + B 2 y l [S+y- l (S-B)B], (68) 

P = — -Kl-x)y-D(_l-x) 1I2 l (69) 

r 

wh ere the last expression i s valid for the idea l gas EoS in Eq. (0, 
see lMignone et all I2005): R yu et alj d2006l) for other options. 

In ECHO the following three inversion methods are imple- 
mented. 

1. The roots of Eqs. d65ll66l l are found simultaneously via a 
two-dimensional Newton technique. This system requires a 
rather accurate initial guess (provided by the quantities found 
at the previous timestep, at the same grid point) and the in- 
version of a 2 x 2 linear system at each iteration. 

2. At each iteration, we derive x = x(y) from Eq. d65l ) and 
then we find the root of /2OO = Fi[x(y),y} — by a one- 
dimensional Newton scheme. This appears to be the most 
straightforward method, since x — x(y) is just a simple alge- 
braic expression, however in the searching process we must 
ensure the condition x < l and sometimes several iterations 
may be required to solve /2OO = 0. 

3. At each iteration, we derive y = y(x) from Eq. ( I66] ) and 
then we find the root of fi(x) = F\ [x, y(x)] = by a one- 
dimensional Newton scheme. This is a variant of the method 
suggested in Paper II and it can only be applied for EoS 
where p is linear in y, as in Eq. j69l . In this case, the root y is 
found either simply as a ratio of two terms, if S ■ B — 0, or as 
the only positive root of the cubic C(y) obtained multiplying 
Eq. d66l ) by y 2 . This may be achieved either analytically or 
numerically via a nested Newton scheme. The existence of 
only one positive root is guaranteed by the following proper- 
ties: C(0) < 0, C'(0) = 0, C(±oo) = +00. 

In the tests presented in Sect. [4] we always use method 3 with 
the nested Newton procedure to find the root of C(y) = nu- 
merically, since it appears to be rather efficient and robust, espe- 
cially when applied to a Newton/bisection hybrid method ensur- 
ing the search of the solution within given boundaries. In cases 
of smooth flows where Eq. < fT0T > replaces the energy equation 
the inversion algorithm is greatly simplified, since sD is the new 
conservative variable, hence the pressure p = sp y depends on x 
alone and we just need to solve the equation f\{x) = 0. 
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3.3. Characteristic speeds in GRMHD 

The spectral properties of the 1-D GRMHD system in Eq. d54l 
are basically the same as for the corresponding system in 
RMHD. Given the structure of the fluxes it is obvious that, for 
example, the eigenvalues of the Jacobian M x will be of the form 



A x = aA' x -B x , 



(70) 



where A' x is the corresponding eigenvalue in special relativistic 
MHD. Thus, in the 3 + 1 approach the gravity terms do not mod- 
ify substantially the hyperbolic structure of the GRMHD equa- 
tions. Full descriptions of the spec tral de composition of the 1-D 
RMHD system in can be found in Anili d!989b . 

Upwind HLL fluxes, described at step 3, just require the cal- 
culation of fast magnetosonic speeds, and this should be accom- 
plished by solving (for each cell and twice for each direction) a 
quartic polynomial, as already described Paper II. However, an 
approximation of these quantities could be also used in Eq. ((55), 
at a price of sli ghtly higher viscosity. In ECHO we follo w the 
strategy by Gamm ie et al.l d2003l) : iLeismann et alj (|2005), who 
realized that, like in classical MHD, an upper bound for fast 
waves is that corresponding to the degenerate case of normal 
propagation k^b^ = 0, where k M = {-to, k x , 0, 0) is the wave 
four-vector. The dispersion relation reduces then to 



(71) 



where the term in square brackets refers to the component of k M 
normal to and 



a 2 — c 2 + c 2 - c 2 c 2 

The sound and Alfven speeds are respectively defined as 

2 _yp 2 _ b 2 



(72) 



(73) 



ph' " ph + b 2, 
where we have introduced the comoving magnetic four-vector 
b" = F ttlv u v = F(v ■ B)n" + B^/T + T(v ■ B)v M , (74) 
and the invariant quantity in Eq. ( |73l is 

b 2 = b^ = B 2 -E 2 = B 2 I T 2 + (v • B) 2 . (75) 

In the degenerate case an analytical expression for the two fast 
magnetosonic characteristic velocities is found by letting A'* = 
CL>lk x in Eq. (TTTT t: 



IX _ (1 - a V ± V fl2 ( 1 - v 2 ) [( 1 - v 2 a 2 )r xx - ( 1 - a 2 )(v x ) 2 ] 



1 — v 2 a 



2 rl 2 



(76) 



and these upper bounds will be then used also for the general, 
non-degenerate case. Note that the above relation, when plugged 
into Eq. ( |70l , correctly reduces to the 3 + 1 GR formu la for the 
hydrodynamical case when B — (Banvuls et al. 1997). 



3.4. Magnetodynamics 

In the pre sent section we summari ze the equations of magneto- 
dynamics ( Komissarov 2002, 2004) and we discuss the few mod- 
ifications implemented in ECHO for the cor responding GRMD 
module. The recipes by iMcKinnevI (l2006al) . which allow one 
to use the same framework of a GRMHD scheme and simply 
neglect the matter contribution, are here followed. In GRMD 
the fluid quantities disappear and the electric field E should re- 
place them as primary variable, together with B. The equations 



to use should be then the two Maxwell equations Eqs. ((6][7), like 
in electrodynamics. However, here we replace Eq. (|6) with the 
electromagnetic momentum-energy conservation law. Thus, by 
setting T 1 ™ 7^ v » T„ in Eqs. and (0 in the limit of neg- 
ligible plasma inertia and thermal contribution, we find 



0. 



(77) 



This force-free situation is actually common to vacuum electro- 
dynamics as well. However, in a highly conducting plasma we 
assume that there is a frame where the electric field vanishes, due 
to the presence of freely moving charges always able to screen 
it efficiently, just like in the GRMHD approximation. This is the 
reason why magnetodynamics is commonly known as degener- 
ate force-free electrodynamics. If the electromagnetic fields are 
decomposed according to the Eulerian observer in the 3 + 1 ap- 
proach of Sect. 12.21 the condition for the existence of a frame 
where the electric field vanishes is replaced by the two invariant 
conditions 



B - E > 0, E ■ B - 0, 



(78) 



which are valid in GRMHD too thanks to ideal Ohm's law 
Eq. ( [391 . If we still indicate with the unit time-like four- 
velocity of this frame, and v is the associated three-velocity de- 
fined in Eq. ( |2TI ). the usual ideal MHD condition is unchanged 
and the two constraints in Eq. ( 1781 are automatically satisfied. 
In order to close the GRMD system, we thus need to express 
this unknown velocity in terms of the electromagnetic quantities 
alone. The required v turns out to be the drift speed of magnetic 
fieldlines 



v = 



ExB 
~B 2 ~' 



(79) 



All the (G)RMHD definitions in Eqs. OH) to ([39) are still valid if 
one neglects matter contribution, in particular 5 = ExB. Notice 
that due to Eqs. ( [391 and ( f79l the three spatial vectors E, B, and 
v are all mutually orthogonal in GRMD. When the three- velocity 
in Eq. ( f79l is used, the equations for GRMHD remain unchanged 
too. However, the continuity equation Eq. ( f30b is now useless, 
while the energy equation Eq. ( [32) is redundant and may be used 
as an additional check. Notice that, in particular, the treatment 
of the metric terms and of their derivatives in the source part 
remains exactly the same as in GRMHD. 

From a computational point of view, the set of GRMD in 
conservative form is easy to treat. The characteristic speeds are 
two Alfven waves and two magnetosonic waves, moving at the 
speed of light. Thus, the expression needed for the simplified 
Riemann solver employed in ECHO (along the x direction) is 
derived from Eqs. ( 170) and d76l by setting a = 1, that is 



A x ± = +a^f^ x7 -B x . 



(80) 



Furthermore, the inversion from conservative to primitive vari- 
ables is also greatly simplified. The magnetic field still enters 
both as a conservative and primitive variable, hence we need to 
derive the drift velocity v for given S and B. The expression em- 
ployed in ECHO is 



1 

B 2 



(S-B) 

B 2 



B 



(81) 



where the second term takes into account the possible numerical 
errors leading to an initial non-vanishing S ■ B. Notice that the 
above formula is equivalent to first derive the electric field as 
E = -S x B/B 2 and then use Eq. d79l . In this way, our code 
preserves the constraint E ■ B = within machine accuracy. 
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4. GRMHD numerical tests 

In order to test our numerical scheme ECHO, several aspects 
need to be checked. First we want to verify that in spite of the 
UCT algorithm, based on staggered representation of the mag- 
netic field components, the overall scheme is able to preserve 
the nominal high order accuracy of the reconstruction and inter- 
polation routines employed. Hence we propose a new test based 
on the propagation of Alfven waves (in flat space-time), which 
are smooth solutions of the equations and thus suitable for such 
kind of problems. However, to better compare ECHO'S perfor- 
mances against other existing GRMHD codes, we will employ 
ECHO at second order in most of the other numerical test prob- 
lems. Thus, even if higher than second order reconstruction al- 
gorithms will be used, in order to sharpen discontinuities and 
reduce numerical diffusion (in particular MP5), all additional 
corrections to achieve an effective higher order of spatial ac- 
curacy will be sometimes disabled and RK2 will be used for 
time stepping in these cases. We will see that the resulting sec- 
ond order scheme (much simpler to be implemented) is a good 
compromise between efficiency, accuracy, and robustness. The 
other numerical tests considered here are: 1-D and 2-D prob- 
lems to check the code shock-capturing properties (a shock tube 
and the cylindrical blast wave); 1-D accretion onto black holes, 
in Schwarzschild and Kerr metrics, to verify ECHO'S high or- 
der properties in curved space-times too; stability of a thick disk 
(with constant angular momentum and with a toroidal magnetic 
field) around a Kerr black hole as a test in 2-D GRMHD. All 
the problems discussed here will involve the presence of sub- 
stantial magnetic fields with plasma beta (the ratio of thermal to 
magnetic pressure) of order of unity or lower. 

If not differently stated, in all our numerical tests we will use 
a Courant number of 0.5, a y-law EoS with y = 4/3, and we will 
solve the equation for the total energy density U. Grid spacing 
will always be constant (though non-uniform grids are permitted 
in ECHO), so the number of points is enough to specify the grid 
in each direction (a single grid point is assigned to the ignorable 
coordinates). 

4.1. Large amplitude CP Alfven wave 

The first test we propose here is a novel one, not previously 
employed in other works on numerical relativistic MHD to our 
knowledge. It involves the propagation of large amplitude circu- 
larly polarized (CP) Alfven waves along a uniform background 
field Bo in a numerical domain, 1 -D or 2-D, with periodic bound- 
ary conditions. Since the propagating wave is an exact solution, 
as we will see below, the test is very useful to check the accuracy 
(both spatial and temporal) and spectral resolution properties of 
a numerical scheme. This is achieved by measuring the errors 
in the solution after one or more periods compared to the ini- 
tial conditions. Such test was first proposed in our Paper II in 
the case of small amplitudes, where the solution was only an 
approximate one. Here we show how to extend the exact solu- 
tion valid in the non-relativistic case to the most general case of 
large amplitudes in (special) relativistic MH D. For the ge neral 
properties o f Alfve nic modes in RMHD see lAnild d 19891) and 
Komissarqv| dl997l) . for other (but less straightforward) numer- 
ical te^tsJnvob/ing_a^ifferejUkind_rf exact solutions 
see iKomissarovl (1 19991) and lDuez etail (120051) . 

Let us consider a CP Alfven wave of normalized ampli- 
tude 77. In classical MHD the variable quantities are the trans- 
verse components of B and v, which are parallel to each other 
with vector tips describing circles in the plane normal to Bq. 
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1-D 




Z-U 




Method 


N 


Li error L 


1 order 


L[ error L 


1 order 


MC2 


8 


1.58e-l 




1 .8 le-1 






16 


3.63e-2 


2.12 


4.60e-2 


1.98 




32 


7.14e-3 


2.34 


8.23e-3 


2.48 




64 


1.55e-3 


2.20 


1.71e-3 


2.27 




128 


3.69e-4 


2.07 


4.01e-4 


2.09 




256 


8.98e-5 


2.04 


9.76e-5 


2.04 




512 


2.21e-5 


2.02 






CEN03 


8 


8.25e-2 




1.07e-l 






16 


1.25e-2 


2.72 


1.68e-2 


2.67 




32 


1.65e-3 


2.92 


2.21e-3 


2.92 




64 


2.09e-4 


2.98 


2.80e-4 


2.98 




128 


2.62e-5 


3.00 


3.50e-5 


3.00 




256 


3.28e-6 


3.00 


4.38e-6 


3.00 




512 


4.10e-7 


3.00 






WEN05 


8 


3.91e-2 




4.76e-2 






16 


2.35e-3 


4.06 


3.14e-3 


3.92 




32 


8.73e-5 


4.75 


1.16e-4 


4.76 




64 


2.82e-6 


4.95 


3.76e-6 


4.95 




128 


8.96e-8 


4.98 


1.19e-7 


4.98 




256 


2.79e-9 


5.01 


3.71e-9 


5.00 




512 


8.53e-ll 


5.03 




_ 


MP5 


8 


1.05e-2 




1.37e-2 






1 ft 
1 u 




A R? 


A QQ P A 


A 78 




32 


1.20e-5 


4.95 


1.16e-5 


4.95 




64 


3.82e-7 


4.97 


5.08e-7 


4.99 




128 


1.20e-8 


4.99 


1.59e-8 


5.00 




256 


3.75e-10 


5.00 


4.98e-10 


5.00 




512 


1.21e-ll 


4.95 







Table 1. Accuracy for the CP Alfven wave test. The L\ errors 
and orders are shown for various methods as a function of the 
number of grid points, both in 1-D and 2-D. Notice that only 
when the error becomes lower than ~ 10" 10 (the value of the tol- 
erance in the inversion from conservative to primitive variables) 
discrepancies from the nominal order start to appear. 



Whatever the wave amplitude, these are the only fluctuating 
fields and the background quantities are not affected by the wave 
(in particular p and p, since the wave is incompressible). In the 
RMHD case, let us look for an exact solution with the same prop- 
erties. The transverse components of B are written 

By = 77B0 cos[k(x - VAt)], B z = 77B0 sin[£(x - v A fj\, (82) 

where we have assumed B x = Bq, va is the (still unknown) 
Alfven speed, and k is the wave vector. Since the induction equa- 
tion remains exactly the same as in the non-relativistic case, we 
still take the velocity components in the form (let us take v x = 
for simplicity): 

v y = -v A B y /B , v z = -v A B z /B , (83) 

as in the classical MHD, where in that case va = Bq/p 1/2 what- 
ever the wave amplitude r\ (the minus sign gives propagation in 
the positive x-direction). We will now see that in the relativis- 
tic case this value is different, basically due to the contribution 
of the kinetic and electromagnetic energies to the inertia of the 
plasma and to the presence of no longer negligible electric forces 
in the momentum equation. 

The electric field is derived from Eq. d39l l, so E y = -v z B x = 
vaB z , E z = v y B x = -VABy, E x = —v y B z + v z B y = 0. Notice also 
that the quantities v 2 = T] 2 v A , B 2 - B 2 } (1 + rj 2 ), and E 2 = rfv 2 A B\ 
are constant, as well as p and p (hence h too). It is easy to show 
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Fig. 1. The large amplitude CP Alfven wave test in the 2-D case (propagation along the diagonal). The L\ errors for the v z velocity 
component, obtained by comparing the solution at the final time t with respect to the initial conditions, for different interpolation 
schemes. In the left panel we show the dependence on the number of grid points N (for a fixed wave number k = 1), whereas in the 
right panel we show the dependence on k for a fixed resolution of 128 x 128. The dashed-dotted line in this second plot refers to the 
run with MP5 at overall second order, and it roughly corresponds to a straight line with L\ ~ k 3 . 



that the transverse components of the momentum equation yield 
the condition 



[ph + {\ 



(84) 



where in square brackets we have the total enthalpy ph + B 2 -E 2 , 
which depends on va itself. Eq. ([84-b is a second order algebraic 
equation for v\, where in order to preserve the condition v 2 A < 1 
the smaller solution must be chosen. Rearranging the terms we 
finally find 



ph+Bld+r, 2 ) 



1+ \ 



2t,B\ 



ph+B 2 (l+T] 2 )) 



(85) 



Notice that in the small amplitude limit 77 <k 1 we retrieve the 
familiar expression v 2 A - B^/(ph + Bh used in Paper II. When 
we further have ft < 1 and fij « p the classical MHD limit 
v 2 A = Bq/p is found, as expected. 

From a numerical point of view, we test the accuracy of our 
scheme by measuring the errors on one of the transverse quanti- 
ties, say v z , at time t = T = L/\>a (one period), compared to the 
initial condition in Eq. ( l82l at t — 0. For the 1-D case we take a 
periodical numerical domain along x of length L - 2n, while in 
the 2-D case we rotate the initial conditions in the (x, y) plane so 
to have propagation along the diagonal of a bi-periodical [0, 27r] 2 
domain. As discussed in Paper II, now two complete spatial peri- 
ods are contained along the diagonal of length L — 2n y2, so we 
can take t = T/2 as final time. With the above choices the wave 
vector k coincides with the wave number, hence it corresponds 
to the (integer) number of spatial periods present in the numeri- 
cal domain. For this test we normalize our physical quantities by 
assuming p = p = Bo = r)=\. 

In Table ([D we show the errors and convergence orders in the 
L\ norm (the absolute error averaged over the whole computa- 
tional domain) for the test with k = 1 at various resolutions. This 
is done for both the 1-D and 2-D cases, and for different recon- 
struction schemes. The errors for the 2-D case are also plotted in 
the left panel of Fig. (fTJ. Note that the nominal order of accuracy 
is achieved already at small numbers N of grid points, which 



means that basically the reconstruction routines employed al- 
ways use the full stencil at their disposal, as expected for smooth 
solutions, without dropping to lower orders at wave extrema. In 
order to achieve third order convergence the RK3 time stepping 
algorithm has been employed for CEN03, while to be able to 
reach an overall fifth order in time and space for WEN05 and 
MP5 the RK3 routine has been used with At oc N~ 5/3 , so that the 
accuracy in time becomes of order (9(Af 3 ) = 0(N~ 5 ) = 0(Ax 5 ), 
i.e. t he same of the spa tial one, as needed in this kind of tests 
(e.g. lJiang & Shulfl996l) . Here the best performing schemes are 
obviously those with higher nominal orders (that for smooth so- 
lutions), thus WEN05 and MP5 in our case. In spite of the same 
fifth order of accuracy, MP5 return smaller errors, up to almost 
a factor 10 at high resolution. This demonstrates that the limit- 
ing conditions in MP5 never apply for this test and the optimal 
stencil is always used, whereas the weights in the WEN05 rou- 
tine do not precisely match to provide the corresponding optimal 
stencil. 

Finally, we test the spectral resolution of our schemes by 
running the same problem at various wave numbers k, from 1 
to 8, at a fixed resolution of 128 x 128 (in the 2-D case). The 
L\ error now increases with k, as expected (an increasing k ba- 
sically means a decreasing resolution), and the dependence is 
stronger for increasing orders r. Indeed, higher order schemes 
(here without the correction to the timestep, thus with a third or- 
der temporal accuracy at most) are able to reproduce reasonably 
the analytical solution even at the smallest wavelengths, where 
second order schemes give poor results. A good compromise be- 
tween efficiency and accuracy is MP5 with RK2 time-stepping 
(3/2 times faster than RK3) and without higher order corrections 
(with overall r — 2 second order accuracy), which appears to be- 
have better than CEN03 with RK3 at small wavelengths. 



4.2. Shock tube with gauge effects 

Shock tubes are excellent tests to monitor the shock-capturing 
properties of a numerical scheme. Until recently, however, an 
exact solver for (special) relativistic MHD Riemann problems 
was still missing, so that comparison was simply made by run- 
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Fig. 2. The relativistic Brio & Wu shock tube modified to allow for gauge effects. The solution on the right hand side refers for a 
run with a — 2, B x = 0, and t = 0.2 (diamonds), whereas that on the left hand side to a run with a — 1, B x — 0.4, and t = 0.16 
(triangles). The numerical solutions are over-plotted to the results obtained with an exact Riemann solver (solid line). Both tests are 
computed with MP5 (no DER and RK2) and N = 1600 grid points. 




Fig. 3. Comparison of different schemes in the relativistic Brio 
& Wu shock tube test. Only N - \QQ grid points are used and 
the density profile is shown for t = 0.4. Results obtained with 
REC based on MP5 (diamonds) appear less smearing than those 
obtained with MC2 (pluses), at the price of some oscillations. 



ning the code at different resolutions and relying on the con- 
vergence properties of the conservative numerical scheme em- 
ployed. Now the situation has changed and we can test our nu- 
m erical solutions against th e exac t Riemann solver for RMHD 
by Giacomazzo & Rezzolla ( 2006), kindly provided by the au- 
thors. Since RMHD shock tubes have been extensively presented 



in Paper II, here just an example will be g iven, namely the rel- 
ativistic version of the Brio & Wu test (lBrio & Wulll988h by 
van Putted d 1993b and lBalsaral (|2001|) . The initial conditions are 



(p,p,B\By) 



( 1.0, 1.0, 0.5, 1.0), x < 0.5 
( 0.125, 0.1, 0.5, -1.0), x > 0.5, 



(86) 



while the other quantities are set to zero. A y-la w EoS with 
y = 2 is used, and the final time is t — 0.4. Following lAnton et al.l 
(2006), instead of showing the standard RMHD results, we turn 
here the test in a sort of GRMHD problem by choosing different 
gauges while preserving a flat metric. In Fig. (O we show the 
numerical results obtained by using a = 2.0 (diamonds), com- 
pared with the exact solution plotted for t/a — 0.2, and those 
obtained with B x = 0.4 (triangles), compared with the exact so- 
lution shifted by 5x = B x t = 0.16. For both runs MP5 is used (no 
DER and RK2), and N = 1600 grid points are employed. 

The first thing to notice is that all the usual structures aris- 
ing from the breakout of the initial discontinuity (left-going fast 
rarefaction wave, left-going slow compound wave, contact dis- 
continuity, right-going slow shock, right-going fast rarefaction 
wave) are well reproduced in both cases, so the chosen gauges 
work as expected. In particular, note the presence near the initial 
discontinuity position x - 0.5 (in the a - 2.0 test) of the so- 
called compound wave, here appearing as a discontinuity. This is 
the combination of an intermediate shock and a rarefaction wave, 
a feature sometimes encountered in coplanar problems due to the 
non-strict hyperbolicity of MHD. Given its nature, it cannot be 
found by exact Riemann solvers and the physical acceptability 
itself as solution of the ideal MHD equations is still de bated 
dBarminet alj|1996t iMvong & Roelll998t iTorrilhonl 120041) . On 
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Fig. 4. The magnetized cylindrical blast wave test at t — 4.0. 2-D maps of density p, thermal pressure p , Lorentz factor F, and 
magnetic pressure p m = B 2 /2 are shown on the four panels on the left hand side. On the two plots on the right we show cuts through 
the center of the domain for thermal (upper panel) and magnetic (lower panel) pressure. Horizonthal cuts are indicated with crosses, 
while pluses are used for vertical cuts. 



the other hand, this feature is invariably found by means of any 
numerical scheme, where some sort of dissipation, either physi- 
cal or numerical, is always present. As far as the reconstruction 
algorithm is concerned, we can see that MP5 gives sharp pro- 
files at all discontinuities, which are captured within 5-10 grid 
points. 

In Fig. (O we show a comparison of the reconstruction 
(RFC) performances of the scheme for the same test, now with 
the original settings (a - l,/3 = and t = 0.4). Here we use low 
resolution runs (N - 100 grid points) to better appreciate the 
differences. The two reconstructions are MP5 and MC2, both at 
an overall second order in space (non DER) and time (RK2). We 
may notice that MP5 provides a more accurate capturing of the 
various waves and discontinuities, in spite of the same overall 
maximum order achieveable, with some extra oscillations, which 
are anyway damped at higher resolutions, as in Fig. [2] Spurious 
oscillations (Gibbs phenomena) near shocks are a well known 
price to pay for high order schemes, especially for those avoid- 
ing decomposition in characteristics, like ECHO. However, we 
deem that the post-processing MP filter behaves quite well in 
this kind of tests. 



4.3. Cylindrical blast wave 

Let us treat RMHD problems involving shocks in more than one 
dimension. A notoriously hard test for relativistic codes is the 
cylindrical blast wave expanding in a plasma with an initially 
uniform magnetic field. This problem was already considered 
in Paper II, here we test our new MP5 scheme a nd we adopt the 
more widely used settings bv lKomissarovl d 1 99 9) . Unfortunately, 
no exact solution is available for the present problem. From a nu- 
merical point of view, in the multidimensional relativistic case it 
is very difficult to treat correctly situations with flow of Alfven 
velocities close to the speed of light. This is because the numer- 
ical errors, which are always present in the reconstruction pro- 
cedures, act independently on, say, x and y components of v and 
B in 2-D runs. This problem easily leads to uncorrect fluxes and 
eventually provides unphysical states, e.g. with v 2 > 1, when 
primitive variables are recovered from the evolved conservative 
ones. Moreover, terms in the total energy equation are strongly 
unbalanced in these cases and, again, numerical errors may lead 
to code crashing. 

The initial conditions are as follows: a square Cartesian box 
[-6, 6] x [-6, 6] contains an internal cylindrical region, within 
r = (x 2 + y 2 ) 1 ^ 2 < 1, with p = 10~ 2 and p = 1. This region is 
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Method 


CPU time 


iter. 


iter. / time 


sub-iter. / time 


MC2-RK2 


86.6 s 


133 


1.54 s~' 


3.08 s" 1 


WEN05-RK2 


97.6 s 


132 


1.35 s- 1 


2.70 s~ l 


MP5-RK2 


100.1 s 


132 


1.32 r 1 


2.64 s- [ 


WEN05-RK3 


152.9 s 


132 


0.86 s- [ 


2.58 a" 1 


MP5-RK3 


154.4 s 


133 


0.86 r 1 


2.58 s- [ 



Table 2. Efficiency results for the cylindrical blast wave prob- 
lem. CPU time (in seconds), total number of iterations, iterations 
per second, and sub-RK cycles per second are reported for vari- 
ous schemes (the DER routine is used only in schemes adopting 
RK3). 



surrounded by an external medium with p = 1CU 4 , p — 5 x 10~ 4 
and these values are reached by means of a smooth ramp func- 
tion betwen r — 0.8 and r = 1. The velocity is zero everywhere 
and the magnetic field is uniform, with B x = 0.1. This is the 
intermediate magnetizat ion case by Komissar ov, with a higher 
external pressure as in iLeismann et al.l d2005l) . We are not able 
to run this test with stronger fields or lower external pressure 
without introducing ad hoc numerical strategies. In Fig. (|4]i we 
show several quantities at t — 4.0, for a run with 200 x 200 
grid points. The scheme used is, as in the previous test, MP5 for 
REC, no DER, and RK2 for time-stepping (overall second or- 
der in both space and time). We notice the presence of several 
structures: an external fast shock, an inner region bounded by 
a reverse shock, both almost circular, and complex anisotropic 
discontinuities in between. Note, in particular, that the magnetic 
field is almost completely swept out from the central region by 
the explosion. The highest outflow speed is reached for y = 
(r max = 3 .69), since there is no magnetic force preventing the ex- 
pansion in the direction along the fieldlines. This problem is also 
a severe test because of the various degeneracies which may oc- 
cur in the Riemann solver. In our case, the HLL procedure with 
the simplified calculation of fast wave speeds does not suffer 
this kind of problems. In spite of the simplified Riemann solver, 
structures appear well defined thanks to the use of an accurate 
REC routine. 

As far as efficiency is concerned, we use the present test to 
measure the CPU time for different scheme settings of ECHO. 
Results are reported in Table (O, where data refer to double pre- 
cision runs on an Intel Xeon 3.0 Ghz processor, for Linux op- 
erating system, with the Intel Fortran compiler. The best per- 
forming scheme is obviously that based on linear reconstruc- 
tion, MC2 in this case, whereas MP5-RK3 is 1.78 times slower. 
However, as we can see from the sub-cycles per second, it is 
the order of the Runge-Kutta method that matters most, whereas 
the DER procedure is quite efficient. When comparing recon- 
struction schemes of the same order, we can notice that MP5 is 
just slightly slower than WEN05 in our implementation, prob- 
ably due to the minmod-type conditions in the limiting process. 
However, from our tests we have found that MP5 is both more 
accurate for smooth solutions and more robust (less oscillatory) 
in problems involving shocks. Our conclusion is that MP5 em- 
ployed at an overall spatial and temporal second order gives the 
best trade-off among efficiency, accuracy and robustness, thus it 
will be used as our base scheme in the next numerical tests. 



4.4. Radial accretion in Schwarzschild metric 

As a first test in a curved space-time we consider here the spher- 
ical transonic accretion onto a non-rotating black hole (of mass 
M = 1) in the presence of a radial magnetic field. The aim is to 
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Fig. 5. Results for the 1-D accretion flows in Schwarzschild 
metric. Quantities are shown by plotting the numerical results 
at t — 100 (diamonds) over the respective exact solution (solid 
line). A resolution corresponding to 100 grid points and recon- 
struction with MP5 are used. 
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Fig. 6. Errors for the 1-D accretion flow in Schwarzschild metric 
of Fig. (0. The L\ norm of the density is shown as a function 
of the grid points N, for MC2-RK2 (diamonds) and MP5-RK3 
(squares). The plasma beta at r c is also varied, according to the 
parameter k — - log 10 /? c . 



check the code ability to preserve in time an analytical solution 
in a curved geometry, where metric terms and their derivatives 
are involved. A full descript i on of t he (fluid) transonic stationary 
solution is given inlMichell d 1972b . here we follow the setup of 
I Anton et"alj 12006). We hence adopt Schwarzschild metric and 
coordinates, with a singular horizon for ry, = 2, where the lapse 
function a = (1 - 2/r) 1 ^ 2 vanishes and y n = a~ 2 diverges. The 
numerical domain is 2.3 < r < 10, the critical point radius 
is r c - 8, an isentropic condition is assumed, and the remain- 
ing free constants are chosen by setting p c = 1/16 (in order to 
have a mass flux of r 2 pTv r = -1) and by assigning the value 
of the plasma beta at the critical radius, B c = 2p c /B 2 , which 
we leave as a free parameter. Note that from an analytical point 
of view the Michel solution does not change in the presence of 
a monopole magnetic field, thus the fluid quantities are unaf- 
fected by the value of the plasma beta (only the magnetic field 
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B' will depend on it, namely as B c ' ), whereas numerically the 
presence of a large magnetic field may lead to severe errors and 
code breaking. This is mainly due to the fact that the numerical 
derivatives of magnetic terms in fluxes do not balance exactly 
the corresponding source terms in the momentum equation, and 
secondly because of the difficulties encountered in the inversion 
routine for the primitive variables. 

In Fig. © we show the results of a simulation with B c - 1 
and N = 100 grid points in the radial direction, comparing the 
quantities obtained at t — 100 with the analytical solutions. 
The scheme employed is MP5 at second order of overall accu- 
racy. Small discrepancies can be seen only near the inner radius, 
where gradients are the largest. To remove both these large gradi- 
ents and the singularity at r — 2 horizon-adapted coord inate sys- 
tems could also be used (Papadopoulos & Font 1998), but here 
we prefer to use the standard Schwarzschild coordinates. For a 
more quantitative comparison, we report in Fig. (|6]l the normal- 
ized L\ errors of the density as a function of the grid points from 
N = 100 to N - 800, for the two schemes MC2-RK2 and MP5- 
RK3 (here with DER). The value of the plasma beta is also var- 
ied, from 1 to 10~ 8 (for an increasing magnetization <x = B 2 /p, 
approximately from 10 _I to 10 7 ). The first thing to notice is that 
the expected scaling with N works also in this non-Cartesian 
case (though there is the usual saturation effect around 10~ 10 ). 
Then we see that the Runge-Kutta order is not an issue in this 
kind of test, where stationary flows are involved (otherwise the 
maximum order would have been 3). Moreover, high resolution 
schemes allow us to reach much lower betas (for N = 800 down 
to B c = 10~ 8 with MP5 at full spatial accuracy order r = 5, and 
B c = 10~ 6 with MC2). If MP5 is employed at second order, inter- 
mediate results are found (not reported in the plot). In order to be 
able to reach such low plasma betas, we have here used Eq. dlOt . 
the adiabatic equation for the entropy function s = p/p r (the so- 
lution is smooth). If the full energy equation is used errors are 
larger by a factor a 2. 

4.5. Equatorial accretion in Kerr metric 

As another example of 1-D test in a curved space-time, we pro- 
ceed further in the level of complexity by studying an accre- 
tion problem in Kerr metric, where not only the lapse func- 
tion a is involved, but also the shift vector B. The problem 
is the magnetiz e d equ atorial flow in Kerr metric described by 
Takahashi et alj d 19901). It is basica lly the general relativistic 
analog of the [Weber & Davis! d 19671) model for the solar wind, 
where the radial velocity has to pass smoothly three critical 
points (slow, fast and Alfvenic) in the equatorial plane where 
the Parker spiral of magnetic field lies. The accretion solution 
was later specialized to the region bet ween the black hole hori- 
zon and the marginally stable orbit bv lGamrni e Hl999T> . in which 
a cold inflow has to cross just the Alfvenic critical point (co- 
incident with the magnetosonic fast point for vanishing thermal 
pr essure). For ou r nume rical test we use the settings prop osed 
bv lGammie eTai] d2003l) and lDe Villiers & Hawlevl l2003), that 
is we study the accretion onto a Kerr black hole with a = 0.5, 
which gives an event horizon at r/, = 1 + (1 - a 2 ) 1 ' 2 1.866 (the 
spherical surface where y„ diverges) and a marginally stable or- 
bit at r mso a; 4.233. After choosing the other free parameters, the 
critical point is located at r c - 3.617. The pressure is initialized 
with an isentropic law, preserving a vanishing thermal contribu- 
tion p <K p => h =i 1 . 

In this test we adopt Boyer-Lindquist coordinates and a ra- 
dial domain 2.1 < r < 4.0 with N = 100 grid points. The results 
are shown in Fig. (0, where the significant physical quantities 
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Fig. 7. Results for the 1-D accretion flow in Kerr metric. 
Quantities are shown by plotting the numerical results at t = 100 
(diamonds) over the respective exact solution (solid line). A res- 
olution corresponding to 100 grid points and reconstruction with 
MP5 are used. 

are plotted at the output time t — 100 against the initial solution. 
As in the previous case, at the outer boundary, where the inflow 
is originated, all quantities are kept constant in time. The scheme 
employed is MP5 at overall second order (no DER and RK2), 
which is rather accurate already at this low resolution, even at 
the inner boundary which is close to the event horizon. For a 
quantitative comparison with the other reconstruction schemes, 
we report here the LI errors on the normalized density as in the 
previous section, again for N - 100 and t = 100. MC2 gives 
2.76e-3, CEN03 (r = 3 and RK3) gives 2.42e-4, both MP5 and 
WEN05 (r = 5 and RK3) give 1.40e-4, while MP5 at second 
order gives 3.22e-4. The improvement of high order methods is 
not as apparent as in the previous test, due to limited precision 
in the initializing routines. 

4.6. Axisymmetric torus in Kerr metric 

The final GRMHD test proposed here is to study the stabil- 
ity of a constant angular momentum thick disk around a Kerr 
black hole and threaded by a toroidal magnetic field. This will 
be achieved through simulations in a 2-D domain, assuming ax- 
isymmetry. For the analytical theory of equili brium in the purely 
hydrodynamical case the reader i s referred tolAbr amowicz et alJ 
d 19781) : iKozlowski et"ail d 19781) : iFont & Daignd d2002l) . while 
the GRMHD versio n with the addition of a p urely toroidal mag- 
netic field is due to Komissar ov et al.l (12006). This test may also 
represent the basis for studying a class of relevant astrophysical 
problems, since the dynamics of accretion disks orbiting around 
black holes is believed to be strongly influenced by the presence 
of magnetic fields. We summarize her e the main features of the 
equilibrium model, while addressing to Komissarov et al. ( 2006) 
for a more detailed description. Under the assumptions of purely 
toroidal velocity and magnetic field, the Bernoulli-like equation 
that needs to be solved is 



dln(-M,) 



Qd£ dp d(R 2 p m ) 



1 - m ph R 2 ph 



= 0, 



(87) 



where I = -u^/u t is the specific angular momentum, £2 = u^/ii' 
is the angular velocity, p m = B 2 /2 is the magnetic pressure 
(notice that the electric field vanishes since v || B), and R 2 = 



L. Del Zanna et al.: ECHO: an Eulerian Conservative High Order scheme for GRMHD and GRMD 



15 



p(x,z) 





0.0 0.2 0.4 



0.8 1.0 






n/4 n/2 3n/4 



0.000 0.005 0.010 0.015 0.020 0.025 



Fig. 8. The results of the magnetized disk evolution. Density (upper panels) and toroidal field (lower panels) are displayed at the 
final output time t = 200. 2-D maps (cylindrical coordinates X = r sin 9, Z - r cos 6 are used for ease of graphical presentation), 
radial and latitudinal cuts through the disk center (r c , tt/2) are shown for the two quantities. Solid lines in the 1-D profiles refer to 
the initial analytical solutions. 



(gt<p) 2 — gttg<M> i s tne generalized distance from the rotation axis. 
We then assume a constant distribution of the specific angular 
momentum, i.e. t = £q, such that Eq. ( T87T ) provides the potential 

R 2 



W = ln(-w t ) = ^ In 



(88) 



The equation of state is barotropic and it is convenient to choose 
p oc (ph) y for the thermal contribution and similarly R 2 p m oc 
(R 2 ph) 7 for the magnetic pressure. Under these assumptions 
Eq. d87l > can be integrated as 



W - W in + 



7 P+Pn 



= o. 



(89) 



y — 1 ph 

The disk is characterized by the condition W < Wi n , where W m is 
calculated at the inner disk radius r m on the equatorial plane. The 
cusp and the center of the disk are defined as those points, again 
in the equatorial plane, where the specific angular momentum 
retains its Keplerian value. Here we use the radius of the disk 
center, r — r c , to determine (q. Notice that the potential W has a 
local minimum at r c , though only in the purely hydrodynamical 
case this point also corresponds to the maxima of p and p. The 
overall disk structure is then completely specified by the two 
radii r m and r c , and by the density p c and plasma beta f3 c at the 
disk center. 

Outside the disk (W > Win) we define a static, unmagnetized 
atmosphere in equilibrium with gravity. This can be obtained by 



adopting the solution of the relativistic Bernoulli equation for an 
isentropic plasma p oc p~* => dp/(ph) = dln/z, and Eq. ( l87l i is 
readily integrated to give the simple relation 



h (-u,) = h (R z /g<p<p) 1/2 = const. 



(90) 



The exact solution can be determined by providing the values 
p atm and p atm , calculated for example at the disk center. Notice 
that all the above relations are valid for both Boyer-Lindquist co - 
ordinates and Kerr-Schild coordinates (e.g. Komi ssaro v 2004), 
which have non-vanishing g tr , g t $ and g rt p terms needed to re- 
move the (unphysical) singularity at the event horizon. In the 
case of Boyer-Lindquist coordinates, employed here for the nu- 
merical test, we have R 2 = a 2 g^ => — u t - a, so that the equilib- 
rium condition for the static atmosphere is simply ah = const. A 
more physical option for the external environment would be to 
define the spherically symmetric Michel's transonic inflow and 
let the system relax to steady state. However, since here we are 
mainly interested in the stability of the disk itself, we prefer to 
use the above static solution, which is an exact one and it is much 
simpler to be initialized. 

The simulation setup is as follows. The numerical domain is 
taken to be 2 < r < 10 and < 6 < n, with 200 grid points in the 
radial direction and 100 in the polar angle direction. We keep the 
quantities fixed in time at both radial boundaries, while reflect- 
ing conditions are imposed at the poles. The first condition is 
needed because otherwise numerical errors near the inner radial 
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boundary, where the gradients are the largest, tend to destabi- 
lize the whole atmosphere. This problem could be also cured by 
choosing appropriate non-uniform grids with higher resolution 
at small radii and/or by using Kerr-Schild coordinates, but here 
we want to retain the simplest possible test conditions. The free 
parameters are chosen to be a = 0.99, r m = 3, r c = 5, p c = 1, 
Pc = 1, Patm = 10~ 5 po p atm = 0.1p atm . Note that the value of 
p c is arbitrary since we are not evolving the metric, which is de- 
termined by the central black hole mass (here taken as unity) 
and angular momentum alone. With the present values we find 
A) 2.80, Wi n - -4.16 x 10~ 2 , W c -9.83 x 10~ 2 so that the 
inner disk disk is located beyond the cusp point and there is a 
finite outer radius (beyond the computational domain). The rota- 
tion period at the disk center is 2w/Q^(r c ) = Inir'J 1 + a) 76.5, 
we take t = 200 as the final output time, corresponding to just a 
few orbital periods but a much longer time with respect to the lo- 
cal dynamical timescales. Here we use MP5 reconstruction at an 
overall r — 2 accuracy (no DER and RK2 for time integration). 

The results are shown in Fig. ([8]), where in the upper pan- 
els we show the density (2-D map, radial cut through the disk 
center, latitudinal cut through the disk center) and in the lower 
panels the toroidal field B*. In the maps we show color images 
and contours of the quantities evolved at the final time (indis- 
tinguishable from those at t — 0), whereas for the 1-D cuts we 
plot the numerical solution at t = 200 (diamonds) together with 
the initial conditions (solid line). Note that our reconstruction 
scheme based on MP5 behaves very well. Minor discrepancies 
appear only near the steep boundaries between the disk and the 
external atmosphere (where density jumps of order 10 3 - 10 4 are 
initially captured by just 2-3 points). Angular momentum is also 
transferred to the non-rotating external atmosphere due to nu- 
merical diffusion in the vicinities of the disk boundaries. The L\ 
norm is 2.8e-4 for the density and 2.41e-5 for the magnetic field, 
while the norm (the largest error in absolute value) is 1 .60e-2 
and 1 .3 le-3, respecti vely. These results a ppear to be comparable 
to those presented by Komissarov (2006a), in spite of the use of 
a much simpler Riemann solver, a constant radial grid spacing, 
and retaining the same overall second order accuracy. Finally, 
errors around the disk boundaries due to numerical diffusion are 
much larger if MC2 is employed instead of MP5, confirming 
that reconstruction based on large stencils may help even near 
discontinuities. The situation is improved if a non-linear radial 
grid is employed, in that case also MC2 provides a good accu- 
racy. On the other hand, results obtained with the full fifth order 
scheme (and RK3 for time-stepping) are similar for this case. 
Finally, note that the present test has been performed by solving 
the full energy equation, and no appreciable changes are noticed 
when Eq. ( fTUb is solved instead. 



5. GRMD numerical tests 

In the present section we perform a series of tests to check the 
performances of ECHO when configured for special and gen- 
eral relativistic magnetodynamics. The numerical settings are 
the same as in the base scheme used for the GRMHD tests, 
namely we employ the HLL solver coupled to MP5 for the re- 
construction (switching off the additional corrections to achieve 
effective higher accuracy) and RK2 for time integration. 

5. 1. Propagation of waves and discontinuities 

Several 1-D tests have been proposed for special relativistic MD. 
Here we select four of them and we change slightly some of the 




Fig. 9. The set of four magnetodynamics 1-D test problems se- 
lected in Sect. 15.11 All plots refer to the B y (x) transverse field 
component at the final output time corresponding to each test. 
From above to below, the four test problems are: fast wave, sta- 
tionary Alfven wave, three-waves, current sheet. 

original setups found in the literature in order to make the no- 
tation more uniform. In all runs we assume a numerical domain 
of 200 grid points in the interval -1.0 < x < 1.0 and a constant 
background field B x - 1 .0. The results of the corresponding sim- 
ulations, at different output times, are all plotted in Fig. [9] where 
the transverse component B y (x) is shown. 

- (a) Fast wave. Her e B z = E x — E y = 0.0 and, following 
iKomissarovl d2002l) . the transverse magnetic field component 
is 



C 1.0, x < -0.6 

B y {x) =\ 1.0 - 1.5(x + 0.6), -0.6 < x < -0.4 
0.7, x > -0.4, 



(91) 



whereas E z (x) = l-B y (x). The fast wave is initially centered 
at x — -0.5 and then should propagate with unchanged pro- 
files at the speed of light. We use t — 1 .0 as output time, so 
that the final position will be x — 0.5. Some small wiggles 
are barely visible in the numerical solution near the corners 
of the wave profile, otherwise the agreement with the analyt- 
ical solution is very good. 

(b ) Stationary Alfven wave. An initial setting similar to that 
in IKomissarovl (|2004|) is assumed for this test, though we 
swap the role of the transverse electromagnetic components 
and the wave profile to make it more similar to the previous 
test. Here we take B z = E y = 1.0, E z = 0.0, and 



[ 1.0, x < -0.1 

B y (x) = \ 1.0 + 1.5(* + 0.1), -0.1 < x < 0.1 
1.3, x > 0.1, 



(92) 



now with E x (x) = -B y (x). This solution is a stationary lin- 
early polarized MD Alfven wave centered at x = 0, so its 
profile should be preserved in time and only numerical dis- 
sipation effects should be found. The output time used for 
this test is t = 2 and from the plot we can see that numerical 
dissipation is negligible for this test, as the initial and final 
profiles are indistinguishable. 

(c) Th ree-waves. This test was proposed by Komissarov 
(2002) and it is concerned with the splitting of a disconti- 
nuity initially located at x = into three waves: two oppo- 
sitely propagating fast waves (traveling at the speed of light) 
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and a standing Alfven wave. It is thus a MD analogue of a 
RMHD shock tube, with the only difference that shocks are 
not allowed in the MD limit. The initial conditions are 



t=0 



t= 1 00 



(B,E) = 



(1.0, 
(1.0, 



1.5, 
2.0, 



3.5, 
2.3, 



-1.0, 
-1.5, 



-0.5, 
1.3, 



0.5), 
-0.5), 



x < 0, 
x> 0, 



(93) 



and the output time is t — 0.75. From the plot we can see that 
the fast wave fronts are reasonably sharp and the the Alfvenic 
discontinuity is preserved within only four grid points. The 
combination of our simple two-waves HLL solver with high 
resolution reconstruction methods like MP5, even when em- 
ployed in an overall second order scheme, confirms thus its 
validity in this kind of tests. Note in particular the absence 
of spurious oscillations, a possible drawback of reconstruc- 
tion methods based on large stencils in problems with sharp 
discontinuities, which proves the limiting capabilities of the 
monotonicity preserving algorithm. 

(d) Current sheet. A current sheet is easily set up by choosing 
B z = E x = £ v = E z = 0.0 and 



B y {x) = 



0.5, 
-0.5, 



x < 0, 
x > 0, 



(94) 



as in [Komissarov (2004). With this value of the transverse 
field, the constraint B 2 - E 2 > is easily preserved through- 
out the evolution. At the output time t = 0.75 we can see 
in the figure the two oppositely propagating fast wave fronts 
located at x = ±0.75, as expected. The numerical diffusion 
of these shocks is very similar to that in the previous test. 

5.2. Uniform magnetic field in Schwarzschild metric 

As a 2-D test in a curved spac e-time we consider the equilibrium 
force-free solution found by Wald (1 19741) . here in Schwarzschild 
metric. An exact solution for the magnetic field is 



B r = B a cos 9, B° = -B ar~ l sin 9, 



(95) 



whereas Z?^ = and E — 0. When translated into cylindrical 
coordinates (R = r sin 9, Z = r cos 9), this is a uniform verti- 
cal field of strength Bq aligned with the Z axis. For our test we 
choose Bq = 1 and a numerical domain 3 < r < 10, < 9 < n, 
with 200 grid points in the radial direction and 100 in the latitu- 
dinal direction. The initial equilibrium is evolved to a large time 
t = 100 (the light crossing time in the radial direction is t = 7) 
and in Fig.[10]we report the magnetic field in vectorial form (the 
length of the arrow is proportional to its strength) for t = and 
t = 100. Only minor discrepancies are visible, for an average 
error of « 6 X 10~ 3 in the field strength. 

6. Conclusions 

We have presented a new code, ECHO, that is the exten- 
sion of our central-type special relativistic scheme (Paper I 
and II) to general relativistic MHD and magnetodynamics. 
This is achieved by applying the general UCT strategies 
dLondrillo & Del Zannafe OOO. 2004} for MHD-like hyperbolic 
systems of conservation laws. The resulting numerical scheme is 
based on simplified Riemann solvers and finite difference high 
order reconstruction methods. As far as the general relativis- 
tic framework is concerned, we adopt here the so-called 3 + 1, 
or Eulerian, formalism. This allows us to present the equations 
(in conservative form) in the most familiar way, i.e. resorting 
to three-dimensional vectors and tensors alone. The limits to 




4 6 
r sin 6 



8 10 




8 10 



Fig. 10. The Wald solution in Schwarzschild metric at the ini- 
tial time t = and output time t = 100. The code employs 
Boyer-Lindquist coordinates, here cylindrical coordinates X = 
r sin 9, Z = r cos 9 are used for ease of graphical presentation. 



special relativistic MHD and classical MHD are then straight- 
forward in this framework. Gravitational terms appear in fluxes 
and in the external sources avoiding the use of complex four- 
dimensional Christoffel symbols. The metric can also depend on 
time and be provided by any solver for Einstein's equations. 

ECHO'S high order procedures are first tested in flat space- 
time, with a new problem involving the propagation of large 
amplitude Alfven waves in 1-D and 2-D domains. We demon- 
strate that the same settings valid for classical MHD can be em- 
ployed in the RMHD problem too by only changing the prop- 
agation speed. This now depends on the amplitude of the wave 
itself, due to the electromagnetic energy contribution to the over- 
all inertia. For the reconstruction routines tested the nominal 
high order of overall accuracy is always reached, up to fifth or- 
der. For the same problem, spectral properties are also checked 
for various schemes by investigating the code behavior at small 
wavelengths, where second order schemes usually fail. Moving 
to discontinuous solutions, one magnetized shock tube is tested 
and even in this case our reconstructions based on larger stencils 
seem to provide sharp profiles also on contact-type discontinu- 
ities, where approximate Riemann solvers usually give poor re- 
sults. In 2-D we study the magnetized blast wave problem, where 
difficulties are known to arise when Cartesian grids are used. We 
find that when the Lorentz factor and/or the magnetization are 
too high, then numerical errors (which are independent along 
each direction) may lead the code to crash. 

In curved space-times, we first study the radial accretion onto 
Schwarzschild black holes in the presence of a monopole mag- 
netic field. High order schemes are able to reproduce the analyt- 
ical solution much better, and this allows us to reach a magne- 
tization as high as 10 7 (for typical values of the other parame- 
ters), while TVD-like second order schemes usually start to fail 
around 10 2 - 10 3 . The expected scaling with the accuracy or- 
der is also reproduced for this test in non-Minkowskian metric. 
In Kerr space-time we test the 1-D equatorial accretion and the 
2-D stability of a constant angular momentum thick disk with 
a toroidal magnetic field (a recently obtained exact solution). 
The latter test provides a very important astrophysical scenario, 
since magnetized tori and rotating black holes are the likely in- 
gredients for AGN and microquasar energy release. Our scheme 
with limited reconstruction based on a five-point stencil is able 
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to maintain the equilibrium solution for several rotation periods 
with negligible errors. This result is achieved without adopting 
specifically designed non-linear grids or horizon adapted coor- 
dinate systems, thus proving its robustness in complex situations 
of astrophysical interest. 

With only minor modifications the GRMHD scheme has 
been also tested in the force-free, lo w-inertia limit of (ideal) 
magnetodynamics (Komissarov 2002). The fluid velocity is re- 
placed by the drift velocity of magnetic field lines and the same 
conservative approach is kept unaltered (McKinnev 2006a). We 
then study the propagation of MD waves and discontinuities 
in flat space-time and the stability of a uniform magnetic field 
around a Schwarzschild black hole. 

As far as efficiency is concerned, a scheme which is fifth 
order accurate in space and third order in time is a: 1.8 times 
slower than a TVD-like second order scheme (in a 2-D test with 
typical resolution). However, most of the difference is due to the 
use of a higher order Runge-Kutta time-stepping algorithm in the 
first case. If the number of sub-cycles per unit time are measured 
instead, then the ratio decreases to just w 1 .2, and therefore high 
order procedures appear to be implemented in a quite efficient 
way. The extra coding necessary to include such routines in an 
existing second order code is not heavy: basically the stencils 
needed for reconstruction must be enlarged, and before every 
derivativation fluxes must be corrected with an additional high 
order (1-D) procedure. 

Concluding, thanks to the Eulerian approach applied to the 
UCT method, we have developed a unified numerical framework 
for MHD-like conservation laws, valid from the classical case 
to special/general relativistic MHD/MD, working in any set of 
curvilinear (even non-orthogonal) coordinates. The base scheme 
conserves the solenoidal constraint for the magnetic field alge- 
braically, due to the UCT strategy, and may be extended to any 
formal accuracy order (for smooth solutions) with finite differ- 
ence upwind reconstruction routines of different kinds. In par- 
ticular, we have proposed her e a limited (filtered) rec onstruction 
based on a fifth order stencil (Sur esh & Huvnhl 1997b . which has 
proved to be both accurate and robust in all the tests performed. 

Acknowledgements. We sincerely thanks Marco Velli, Simone Landi, Antonio 
Scalabrella for fruitful discussions, and an anonymous referee for the precious 
suggestions which have helped us to improve the scientific quality of this work. 
Olindo Zanotti was supported by MIUR through COFIN funds (Pacini). Niccolo' 
Bucciantini was supported by NASA through Hubble Fellowship Grant HST- 
HF-01 193.01A, awarded by the Space Telescope Science Institute, which is op- 
erated by the Association of Universities for Research in Astronomy, Inc., for 
NASA, under contract NAS 5-26555. 



References 

Abramowicz, M., Jaroszynski, M., & Sikora, M. 1978, A&A, 63, 221 
Aloy, M. A., Ibanez, J. M., Marti, J. M., & Muller, E. 1999, ApJS, 122, 151 
Anile, A. M. 1989, Relativistic fluids and magneto-fluids (Cambdridge 

University Press) 
Anninos, P., Fragile, P. C, & Salmonson, J. D. 2005, ApJ, 635, 723 
Anton, L., Zanotti, O., Miralles, J. A., et al. 2006, ApJ, 637, 296 
Arnowitt, R., Deser, S., & Misner, C. W. 1962, The Dynamics of General 

Relativity, ed. L. Witten (New York: John Wiley), 227-265 
Balsara, D. 2001, ApJS, 132, 83 

Balsara, D. S. & Shu, C.-W. 2000, Journal of Computational Physics, 160, 405 
Banyuls, F, Font, J. A., Ibanez, J. M., Marti', J. M., & Miralles, J. A. 1997, ApJ, 
476, 221 

Barmin, A. A., Kulikovskiy, A. G., & Pogorelov, N. V. 1996, J. Comput. Phys., 
126,77 

Baumgarte, T. W. & Shapiro, S. L. 2003, ApJ, 585, 921 
Blandford, R. D. & Znajek, R. L. 1977, MNRAS, 179, 433 
Brio, M. & Wu, C. C. 1988, J. Comput. Phys., 75, 400 
Bucciantini, N., Amato, E., & Del Zanna, L. 2005, A&A, 434, 189 



Bucciantini, N., Thompson, T. A., Arons, J., Quataert, E., & Del Zanna, L. 2006, 

MNRAS, 368, 1717 
Colella, P. & Woodward, P. R. 1984, Journal of Computational Physics, 54, 174 
De Villiers, J.-P. & Hawley, J. F. 2003, ApJ, 589, 458 
De Villiers, J.-P., Hawley, J. F, & Krolik, J. H. 2003, ApJ, 599, 1238 
De Villiers, J.-P., Hawley, J. E, Krolik, J. H., & Hirose, S. 2005, ApJ, 620, 878 
Del Zanna, L., Amato, E., & Bucciantini, N. 2004, A&A, 421, 1063 
Del Zanna, L. & Bucciantini, N. 2002, A&A, 390, 1 177 
Del Zanna, L., Bucciantini, N., & Londrillo, P. 2003, A&A, 400, 397 
Del Zanna, L., Volpi, D., Amato, E., & Bucciantini, N. 2006, A&A, 453, 621 
Duez, M. D., Liu, Y. T, Shapiro, S. L., Shibata, M., & Stephens, B. C. 2006a, 

Physical Review Letters, 96, 031101 
Duez, M. D., Liu, Y. T., Shapiro, S. L., Shibata, M., & Stephens, B. C. 2006b, 

Phys. Rev. D, 73, 104015 
Duez, M. D., Liu, Y. T., Shapiro, S. L., & Stephens, B. C. 2005, Phys. Rev. D, 

72, 024028 

Eulderink, F. & Mellema, G. 1994, A&A, 284, 654 
Evans, C. R. & Hawley, J. F. 1988, ApJ, 332, 659 
Font, J. A. 2003, Living Reviews in Relativity, 6, 4 
Font, J. A. & Daigne, F. 2002, MNRAS, 334, 383 

Font, J. A., Ibanez, J. M., Marquina, A., & Marti, J. M. 1994, A&A, 282, 304 
Gammie, C. F. 1999, ApJ, 522, L57 

Gammie, C. F, McKinney, J. C, & Toth, G. 2003, ApJ, 589, 444 
Gammie, C. F, Shapiro, S. L., & McKinney, J. C. 2004, ApJ, 602, 312 
Giacomazzo, B. & Rezzolla, L. 2006, J. Fluid. Mech., 562, 223 
Goldreich, P. & Julian, W. H. 1969, ApJ, 157, 869 

Hai1en, A., Engquist, B., Osher, S., & Chakravarthy, S. 1987, J. Comput. Phys., 
71,231 

Hai-ten, A., Lax, P. D., & Van Leer, B. 1983, SIAM Rev., 5, 1 
Hawley, J. F. & Krolik, J. H. 2006, ApJ, 641, 103 
Jiang, G.-S. & Shu, C.-W. 1996, J. Comput. Phys., 126, 202 
Koide, S. 2003, Phys. Rev. D, 67, 104010 

Koide, S., Kudoh, T, & Shibata, K. 2006, Phys. Rev. D, 74, 044005 

Koide, S., Meier, D. L., Shibata, K, & Kudoh, T. 2000, ApJ, 536, 668 

Koide, S., Shibata, K, & Kudoh, T. 1999, ApJ, 522, 727 

Komissarov, S. S. 1997, Physics Letters A, 232, 435 

Komissarov, S. S. 1999, MNRAS, 303, 343 

Komissarov, S. S. 2001, MNRAS, 326, L41 

Komissarov, S. S. 2002, MNRAS, 336, 759 

Komissarov, S. S. 2004, MNRAS, 350, 427 

Komissarov, S. S. 2005, MNRAS, 359, 801 

Komissarov, S. S. 2006a, MNRAS, 368, 993 

Komissarov, S. S. 2006b, MNRAS, 367, 19 

Komissarov, S. S., Barkov, M., & Lyutikov, M. 2006, MNRAS, 1280 
Komissarov, S. S. & Lyubarsky, Y. E. 2004, MNRAS, 349, 779 
Kozlowski, M., Jaroszynski, M., & Abramowicz, M. A. 1978, A&A, 63, 209 
Landau, L. D. & Lifshitz, E. M. 1962, The classical theory of fields (Oxford: 
Pergamon) 

Leismann, T, Anton, L., Aloy, M. A., et al. 2005, A&A, 436, 503 
Lele, S. K. 1992, J. Comput. Phys., 103, 16 

Lichnerowicz, A. 1967, Relativistic hydrodynamics and magnetohydrodynamics 

(New York: Benjamin) 
Liu, X.-D. & Osher, S. 1998, J. Comput. Phys., 141, 1 
Londrillo, P. & Del Zanna, L. 2000, ApJ, 530, 508 
Londrillo, P. & Del Zanna, L. 2004, J. Comput. Phys., 195, 17 
Martf, J. M. & Muller, E. 1996, J. Comput. Phys., 123, 1 
Marti, J. M. & Muller, E. 2003, Living Reviews in Relativity, 6, 7 
McKinney, J. C. 2005, ApJ, 630, L5 
McKinney, J. C. 2006a, MNRAS, 367, 1797 
McKinney, J. C. 2006b, MNRAS, 368, 1561 
McKinney, J. C. 2006c, MNRAS, 368, L30 
McKinney, J. C. & Gammie, C. F. 2004, ApJ, 61 1, 977 
Michel, F. C. 1972, Ap&SS, 15, 153 
Mignone, A. & Bodo, G. 2006, MNRAS, 368, 1040 
Mignone, A., Plewa, T., & Bodo, G. 2005, ApJS, 160, 199 
Misner, C. W., Thorne, K. S., & Wheeler, J. A. 1973, Gravitation (San Francisco: 

Freeman) 

Mizuno, Y, Yamada, S., Koide, S., & Shibata, K. 2004, ApJ, 615, 389 
Myong, R. S. & Roe, P. L. 1998, J. Comput. Phys., 147, 545 
Nishikawa, K.-I., Richardson, G., Koide, S., et al. 2005, ApJ, 625, 60 
Noble, S. C, Gammie, C. E, McKinney, J. C, & Del Zanna, L. 2006, ApJ, 641, 
626 

Papadopoulos, P. & Font, J. A. 1998, Phys. Rev. D, 58, 024005 
Pons, J. A., Font, J. A., Ibanez, J. M., Marti, J. M., & Miralles, J. A. 1998, A&A, 
339, 638 

Roe, P. L. 1981, J. Comput. Phys., 43, 357 

Ryu, D., Chattopadhyay, I., & Choi, E. 2006, ApJS, 166, 410 



L. Del Zanna et al.: ECHO: an Eulerian Conservative High Order scheme for GRMHD and GRMD 



Shibata, M., Duez, M. D., Liu, Y. X, Shapiro, S. L., & Stephens, B. C. 2006, 

Physical Review Letters, 96, 031 102 
Shibata, M. & Sekiguchi, Y.-I. 2005, Phys. Rev. D, 72, 044014 
Shu, C.-W. 1997, NASA ICASE Rep., 97, 65 
Shu, C.-W. & Osher, S. 1988, J. Comput. Phys., 77, 439 

Sloan, J. & Smarr, L. L. 1985, General relativistic magnetohydrodynamics, ed. 

J. L. J. Centrella & R. Bowers (Boston: Jones and Bartlett), 52-68 
Smarr, L. & York, Jr., J. W. 1978, Phys. Rev. D, 17, 2529 
Spitkovsky, A. 2006, ApJ, 648, L51 

Suresh, A. & Huynh, H. T. 1997. J. Comput. Phys., 136, 83 

Takahashi, M., Nitta, S., Tatematsu, Y, & Tomimatsu, A. 1990, ApJ, 363, 206 

Thorne, K. S. & MacDonald, D. 1982, MNRAS, 198, 339 

Torrilhon, M. 2004, Journal of Plasma Physics, 69, 253 

van Putten, M. H. P. M. 1993, J. Comput. Phys., 99, 341 

Wald, R. M. 1974, Phys. Rev. D, 10, 1680 

Weber, E. J. & Davis, L. J. 1967, ApJ, 148, 217 

Weinberg, S. 1972, Gravitation and Cosmology (New York: Wiley) 

Wilson, J. R. & Mathews, G. J. 2003, Relativistic Numerical Hydrodynamics 

(Cambridge, UK: Cambridge University Press) 
Yee, K. 1966, IEEE Trans. Antennas Propag., 14, 302 

York, Jr., J. W. 1979, Kinematics and dynamics of general relativity, ed. L. Smarr 

(Cambridge: Cambridge Univ. Press), 83—126 
Zhang, X.-H. 1989, Phys. Rev. D, 39, 2933 



Appendix A: Finite difference procedures 

ECHO employs finite difference piecewise polynomial high or- 
der procedures for interpolation, reconstruction, and derivation. 
Compact-like (implicit) routines are also implemented, but we 
do not discuss them here. Below we will indicate with n the or- 
der of accuracy of the single procedures, while r will retain the 
meaning of the spatial accuracy of the overall scheme. 

A. 1. Interpolation (INT) 

Interpolation is explicitly needed to approximate the magnetic 
field components at step 1 in Sect. 13.11 but it also provide the 
building blocks for upwind reconstruction methods. For any kind 
of polynomial interpolation, it is convenient to calculate the co- 
efficients by means of the Lagrange formula. For a stencil of n 
points Xi (either cell centers or intercell points), the polynomial 
approximating a function f(x) to n-th order is 

n n 

Pn(x)=/aifi, di= |f , (A.l) 

H k=lMi Xi ~ Xk 

where by construction p n (xf) = /; = f(xf). For the case of mag- 
netic field interpolation, we need to approximate a function f(x) 
at cell center xj for given intercell values /j+i/2- Application of 
Eq. dA.ll ) to a symmetric stencil around xj gives the expressions 

fj = {fj-H2+f j+m )l2, (A.2) 
/; = (-/;-3/2 + 9/^-i/2 + 9f j+l/2 - /; +3 / 2 )/16, (A.3) 

/.. = 0//-5/2 - 25/;_ 3/2 + 150/)- 1/2 + 

+ 150/) +1/2 - 25/) +3/2 + 3/- +5 /2)/256, (A.4) 

respectively for n = 2, n = 4, and n - 6. Thus, the «-th order 
formula should be used for an overall scheme with r <n. 

A.2. Reconstruction (REC) 

The reconstruction process employed in ECHO is again an op- 
eration based on piecewise polynomial interpolation. Given a 
stencil of n grid points [xj] (cell centers) with corresponding val- 
ues {/)} of the discretized function fix) (in ECHO the primitive 
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variables, see step 2), the problem is to find a n-th order ap- 
proximation of the intercell value /}+i/2- Note that in the numer- 
ical literature high order reconstruction is usually implemented 
to find directly the fj+1/2 numerical fluxes of step 4 (called re- 
construction via the primitive function), corresponding to our 
REC+DER combined operations. Therefore, the polynomial co- 
efficients presented here will differ with those usually found in 
the literature. Contrary to the centered interpolations seen above, 
in shock-capturing schemes upwind interpolation is needed, that 
is based on either left-biased (L) or right-biased (R) stencils. To 
achieve this, typically n is chosen an odd number and the two 
stencils are taken symmetric with respect to Xj+i/2- Moreover, 
the same reconstruction routine H({fj}) may be employed for 
both L and R procedures: 

/f+l/2 = W;-(»-l)/2>--->/;+(«-l)/2)> (A.5) 
ff+l/2 = K(//+1+(h-1)/2> • • • . fj+l-(n-l)/2)- (A.6) 

For n — 1 we have the expected upwind constant approxima- 
tions fji l/2 = fj, ff +l/2 = f j+ i. For n > 1 we have to face the 
problem that the two stencils may contain a discontinuity, hence 
sub-stencils should be used in order to avoid Gibbs oscillations 
and the above formula actually refers to the optimal stencils pro- 
viding an order n only for smooth solutions. 

For n — 3 we have quadratic interpolation. By applying 
Eq. ( IA.1I ). the left fixed-stencil reconstruction (only left recon- 
structions will be considered hereafter) based on the optimal 
stencil is 

/>i/2 = i-fj-i + 6fj + 3/,- +I )/8. (A.7) 

In TVD-like reconstructions, based on the same n — 3 sten- 
cil used above, third order is sacrificed for sake of stability by 
resorting to second order for continuos fields and to first or- 
der when a discontinuity is present. These schemes are based 
on piecewise linear reconstruction and monotonicity is typically 
enforced by making use of slope limiters 

//'+1/2 = fj + \S (A_fj,A+fj), (A.8) 

where A ± f> = +{fj±\ - fj) and the slope S can be for example 
the MinMod (MM2 in ECHO) limiter 

mm(x,y) = {-[sgn(x) + sgn(y)]min(|x|, |y|), (A.9) 

or the so-called Monotonized Centered (MC2 in ECHO) limiter 

mc(x,y) = i[sgn(jc) + sgn(v)]min(2|x|,2|y|,i|x+y|). (A.10) 

Usually reconstruction based on MM2 is safer but more smear- 
ing, while MC2 provides a good compromise between robust- 
ness and accuracy. Note that at local (smooth) extrema all lim- 
ited reconstructions of this kind drop to first order. ENO schemes 
follow a different strategy. In EN02, one between the two linear 
interpolations based on 2-point sub-stencils 

fjll/2 = fj + L 2 A -fi - (-//"> + 3/)) A ( A -ll) 

f?V2 = fi + l A +fj - (fj + /*i)/ 2 . ( A - 12 > 

is chosen, with selection procedures based on smoothness cri- 
teria to ensure the (essentially) non-oscillatory behavior. Thus, 
EN02 always employs a piecewise linear interpolation. The 
possibility to achieve the optimal third order reconstruction of 
Eq. iA.71 is provided by the weighting process in the WEN03 
procedure: 

/; + i/2 = ^i// + i I ) /2 +w2/; + 2 1 ) /2 , (a. 13) 
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where the optimal reconstruction is found for a>i = 1/4 and a>2 - 
3/4. In the (nonlinear) selection process these are the limits for 
smooth fields, otherwise a different combination (resulting in a 
lower order) is achieved and for discontinuous fields WEN03 is 
equivalent to EN02. 

Analogous possibilities for ENO-like schemes are offered by 
reconstruction based on the n — 5 stencil. The optimal choice 
yielding fifth order accuracy is 

W = Ofj-i - 20/m + 90// + 60//H - 5/; +2 )/128, (A. 14) 

while the three 3-point sub-stencils provide the quadratic inter- 
polations 

film = ( 3 /j'-2 - 10/M + 15//)/8, (A.15) 
/®/2 = ("//-I + 6 /; + 3/^0/8, (A.16) 
/$/2 = + 6 /i+i - //+2)/8, (A. 17) 

which are easily obtained as usual by either use of Eq. ( IA. lb 
or by Taylor expansion. In EN03 third order reconstruction is 
always obtained by choosing the smoothest among the above 
interpolations. In WEN05 the combination 

W = mf$i2 + ^mp. + <*f?w> (A.18) 

is used, and the optimal fifth order reconstruction in Eq. ( IA.14b is 
retrieved when u)\ = 1/16, a>2 = 10/16 and W3 = 5/16, obtained 
for smooth fields. Another pos sibility is provided by the CEN03 
algorithm dLiu & O sher 1998), which is basically equivalent to 
EN03 for smooth fields (thus both achieve third order at most) 
and it reduces to lower order TVD reconstruction (hence even to 
first order) in the presence of discontinuities, but not at smooth 
extrema. The robustness and accuracy of this scheme were com- 
prehesively tested in Paper I and II. 

A different strategy is followed by M P {Monotonicity 
Preserving) methods dSuresh & Huvnhll 19971) : first the high or- 
der reconstruction, like that in Eq. ( IA. 141 > for MP5, is con- 
structed, then, if spurious oscillations are found, a nonlin- 
ear filter based on limiting algorithms is applied to reduce 
them, retrieving first order approximations only where needed 
(like in CENO). An approach similar to MP is that fol- 
lowed in the celebrated P PM (Piecewise Parabolic Method, 
IColella & Woodward! fl98l) . very popular among astrophysi- 
cists, due to the rather sharp profiles provided at disconti- 
nuities, and used in s pecial relativist i c HP and MHD too 
dMartf & MiiHerl 1 19961: iMignone etail 120051: iLeismann et alj 
2005). However, that method has the drawback of reducing to 
first order even at smooth extrema, just like TVD. Moreover, 
the post-processing filters for PPM are rather involved and heav- 
ily system-dependent(especially the steepening of contact-like 
discontinuities), thus in conflict with the philosophy adopted 
here. On the other hand, MP methods are particularly suitable 
for component-wise reconstruction and these filters can be ap- 
plied to a variety of explicit in terpolants, to higher order WENO 
methods (Bals ara & Shdl2O00r) . or even to co mpact interpola- 
tions with spectral-like resolution dLeldl 1992b . The MP5 algo- 
rithm based on the n — 5 explicit reconstruction of Eq. jA. 14b 
has been shown here to be both highly accurate and robust in all 
tests, and we thus recommend its use. We refer to the original 
paper for a description of the nonlinear filter. 

A.3. Derivation (DER) 

The derivation operation was encountered at step 4 to provide the 
numerical flux function fj+1/2, given a stencil of intercell fluxes 



{//+1/2L This must be done in such a way that (/ J+ 1/2-/7- 1/2) A i s 
an appropriate high order approximation of the fix) first deriva- 
tive calculated at x = Xj, where h is the (constant) grid spacing. 
Let us then start by looking for a finite difference approximation 
of the first derivative. It is convenient to write it as 

hf'(Xj) w $+1/2 - /j-l/2 = fl(/y+l/2 - fj-1/2) + 

+b(fj + 3/2 ~ /y-3/2) + c(// +5 /2 - /j-5/2), (A. 19) 

where we have truncated the approximation up to sixth order. If 
we now expand both sides of the above equation in Taylor series 
around Xj we find the system 




(A.20) 



where the exponents indicate derivation of the corresponding or- 
der and where clearly all terms with even k vanish. For n = 2, 
where b — c — 0, we simply find a — 1 . For n = 4, where only 
c — 0, the above system is readily solved by a = 9/8,/? = -1 /24. 
Finally, for n — 6 the solution is a — 75/64, b = -25/384, 
c = 3/640. The next step is to write 

fj+1/2 = do f j+1/2 + <Wj-l/2 + /y+3/2) + d 4 (fj-3/2 + fj+5/2% (A.21) 

and comparison with Eq. ( IA.19I ) provides the relations do — a + 
b + c, d\ — b + c, ds, — c. For n — 2 do — 1, di = di = 
and fj+1/2 = fj+1/2, as expected. Thus, no extra high order 
corrections on numerical fluxes are needed for schemes up to 
second order. For n — 4 we find do = 13/12, <i 2 = -1/24, and 
d 4 = 0. Finally, for n = 6wefind£/o = 1067/960, d 2 = -29/480, 
and d A = 3/640. 

In order to highlight the nature of the DER procedure as a 
correction for higher than second order approximations, it is con- 
venient to rewrite Eq. dA.21| ) in the form 

fj + H2 = f J+ i/2 ~ ^A% 1/2 + JLAWf j+m , (A.22) 

where only the first term is retained for n = 2, the second is 
introduced for n = 4, and the complete expression is used for n = 
6. For a generic index i the second and fourth order numerical 
derivatives are respectively given by 

A C2) /i = /-1-2/ + /+1, (A.23) 
A< 4 >/ = A( 2 VU-2A®/i+A (2) / + i = 

= Z-2 - 4/ w + 6/ - 4/+i + f i+2 . (A.24) 

Notice that here only DER operators based on centered, sym- 
metric stencils have been considered. The high order corrections 
described above can be easily turned into non-oscillatory algo- 
rithms by any sort of limiting or stencil selection upwind pro- 
cess, like those employed for REC. 



